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ABSTRACT 


The requirement to develop a space based remote ocean 
sensing platform exists within the Department of the Navy. 
This project models a satellite subsystem with structural 
flexibility utilizing the Equivalent Rigid Link System 
(ERLS). Dynamic analysis with computer simulation is 
presented for a simple flexible boom rotating in three 


dimensions with and without a point mass at the boom tip. 
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TABLE I DEFINITION OF TERMS 


KE = Kinetic Energy 

PE = Potential Energy 

8 = Large motion generalized coordinate 

F = Large motion generalized force, applied moment 

U = column vector for small motion generalized 
coordinates. 

= = local position vector (L subscript refers to load) 

R = global position vector (L subscript refers to load) 

qd = deformation vector of transverse deflections 

Wl = mass density of the link (L subscript refers to load) 

W = transformation matrix (a function of 6) : 

Wo = derivative of the transformation matrix wrt 8 

S = residual accelerations term 

D = deformation transformation matrix 

® = shape function matrix derived in Appendix C 

r = spacial derivative of the shape function matrix 

C = stiffness matrix 

q = gravitational acceleration vector 

Loy = volume integral with a kernel for dependent matrices 
where ## represents an index number from Table A 

K, = length integral (along the boom) 

Hy = volume integral where the # represents an index 


number # from Table A.2 
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A. BACKGROUND 

The U.S. Navy has declared it's intention to enhance its 
global, real-time monitoring of the ocean characteristics and 
to improve its ability to forecast oceanographic parameters 
by means of a dedicated space system [Ref. 1]. The system 
under consideration would improve the battle group 
commander's ability to screen his forces from detection by 
enemy submarines, be more aware of severe weather conditions 
in the operating theater and improve forecasters meteorology 
predictions. 

The satellite will use a variety of special sensors and 
Operate from an 820 km polar orbit. A number of physical 
considerations for this type of mission embody the foreground 
of modern structural and control engineering research. 

One such item considered necessary to achieve the Navy's 
desired goals is a Low Frequency Microwave Radiometer (LFMR). 
It represents a unique aspect of spacecraft engineering and 
GOneErOo!. The LFMR would provide sea surface temperatures 
accurate to within a half of a degree Kelvin with a surface 
resolution capability of 25 km. The LFMR would require a 
large (approximately 5 meter diameter) antenna with precise 
pointing knowledge; its field of view would scan at either 31 


or 16 scans per minute. The width of scan would be at 
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approximately 1400 km centered on the ground track, in front 
of or behind the satellite. The incident beam would make a 
53.1 degree angle with the local vertical at the surface, see 
Figure 1. (Ref. 2:p. 3] 

The large ground coverage is provided by the 
eccentrically mounted antenna. A schematic drawing of the 
proposed LFMR sub-system required to support the antenna is 
presented in Figure 2. The boom sizes and off-axis angles 
vary from contractor to contractor, those represented in 


Figure 2 are not considered the only feasible choices. 


B. MOTIVATION 

Currently the Navy .has halted development on the Navy 
Remote Ocean Sensing System (N-ROSS), which was to be the 
vehicle’ for the LFMR. Cost considerations and_ the 
technological development (simulations and modeling ) of the 
LFMR boom assembly and control scheme have been suggested as 
primary concerns. 

This investigation is motivated by the limited number of 
techniques currently employed for the dynamic analysis of 
large flexible space structures. We will present’ the 
development and computer simulation of an analysis for a 


flexible spacecraft boom, rotating in three dimensions. 


is 


C. LITERATURE REVIEW 

The approach to modeling the kinematics and dynamics of 
the rotating flexible space antenna are not unlike those 
techniques employed in the field of robotics necessary for 
the design and control of flexible manipulators. In the case 
of flexible manipulator models, the focus is to develope the 
relationship between the large, rigid-body motion (large 
motion) and the small motions due to the structanea 
flexibility. 

The most popular class of these reduction methods for 
dynamic analysis of structures appears to employ the modal 
Superposition technique [{Refs. 3-7]. In this type of 
analysis a determination must paunnae to limit the number of 
contributing mode shapes and then combine those that will 
contribute in away which yields a system of equations of 
motion. 

Book [{Refs. 4,5] utilizes this modal technique to model a 
system with flexible kinematics. He utilizes recursive 
Lagrangian kinematics and large and small motion coordinates. 
The resulting equations of motion are an accurate model 
describing large and small motion and the coupling between 
them. However the result, a set of coupled, non-linear 
second order, ordinary differential equations, proved 
computationally intensive, and therefore expensive to solve. 

Sunada and Dubowsky [Ref. 8] on the other hand expressed 


the kinematics and dynamics of the manipulator in terms of 
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"Small perturbations about some Known nominal motion of the 
system". The effect of the coupling between the large and 
small motion is incomplete when the large motion is assumed 
nominal and the small motion superimposed upon it. 

Chang's development of the Equivalent Rigid Link System 
ERLS [Ref. 9], describes the large motion kinematics of the 
system by the ERLS and the small motion deformations relative 
Eo iC. The application of finite element techniques and 
Lagrangian dynamics, produces two sets of coupled, non- 
linear, ordinary differential equations of motion. The 
application of the ERLS formulates these sets of equations; 
one set for large motions and one set for small motions. The 
set of large motion equations are non-linear in both the 
large and small motion variables. The set of small motion 
equations are linear in the small motion variable and non- 
linear in the large motion variables. 

The characteristics of these equation sets permit a 
solution by numerical integration with a variety of 
techniques possible. The integration technique developed by 
Chang, the Sequential Integration Method, is straightforward 
and relatively easy to implement. Chang's model offered a 
complete representation of the dynamics of flexible 
Manipulators as well as one that could be efficiently solved 
using the Sequential Integration Method. 

Simo and Vu-Quoc [Ref£. 10] present an approach which 


appears well suited for design but the choice of generalized 
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coordinates may require a complicated observer to relate the 
computational prediction of the generalized coordinates with 
observable reference quantities of the controller 
environment. 

Previous work at the Naval Postgraduate School utilizing 
the ERLS technique validated the procedure in a planar 
application PRetsee) 117121. Both Petroka and Gannon applied 
Chang's ERLS utilizing separate mode shape response matrices. 
The agreement between the assumed modal responses and the 
experimentally observed response justifies this first effort 
at applying the ERLS to model the dynamic response of 


flexible booms in the spacecraft applications. 


D. OUTLINE 

The nonlinear dynamic response of a spacecraft boom 
structure member is predicted by the numerical integration of 
the equations of motion. In this work, the spatial finite 
element discretization of the boom and the application of an 
assumed polynomial mode shape provides an approximate 
solution. 

The modeling process has been applied through the ERLS 
utilizing Lagrangian Dynamics to obtain the system equations 
of motion for numerical integration. 

Chapter Two of this work reviews Chang's development of 
the ERLS and in Chapter Three the ERLS is applied toa 


rotating flexible boom and a variety o£ simulation results 
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are presented. The computer simulation code and methodology 
are introduced in Chapter Four. The results of this work 
illustrate both the need for controlling the flexibility and 
the versatility o£ the ERLS solution technique. A review of 
this work and suggested areas for development is presented in 


Chapter Five. 
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Il. OYNAMIC MODEL FORMULATION FOR A SPACECRAFT 
GLE- M ri IBLE BOOM 


A. KINEMATIC RELATIONSHIPS OF THE EQUIVALENT RIGID LINK 

SYSTEM 

The introduction and development of the Equivalent Rigid 
Link System, separates the Kinematics of the flexible body 
into large and small motions. The large motions are modeled 
with the characteristics of rigid bodies while the small 
motions are modeled by finite element techniques and defined 
relative to the position of the equivalent rigid bodies. 

The small motions resulting from loading flexure are 
- measured in the ERLS relative to the local coordinate systein 
for each equivalent rigid body. The absolute position of a 
point on a deformed body is obtained by coordinate 
transformations which utilize the characteristics of joint. 
variables as in the methods of Hartenberg and Denavit (Ref. 
13). Lagrangian dynamics are used to obtain the equations of 
motion and therefore time derivatives of the absolute 


position will be necessary in the development of the 


Lagrangian expression. 

The ERLS model accounts for the change in location 
experienced by each point along a flexible member. Discrete 
points or nodes are chosen and their relative position 
determined while locations of intermediate points are 


assumed according to a mathematical function. The choice of 
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nodes and mathematical function determines the theoretical 
shape of the flexed member. 

The Finite Element Method (FEM) is used to descretize the 
flexible member assigning the nodes and enforcing 
compatibility between elements. The choice of shape function 
and node position must be consistent with the accuracy 
desired. The ERLS will accommodate closed form natural mode 
shape functions or any order of polynomial. In this work for 
example a cubic polynomial shape function was chosen to 
represent the bending displacement of the flexible boom. The 
nodal positions chosen here are _ the end points of the 
flexible member. The selection of a- polynomial shape 
function and end point nodal positions are consistent with 
previous works which have experimentally validated their 
selection. (Refs. 9,11,12]. 

The kinematics of the ERLS separates the global motions 
into the sums of large and small motions. The large motion 
is described by the equivalent rigid link and the small 
motion due to deformations of the link are described relative 
to the rigid link. 

The characteristics of rigid body dynamics is essentially 
classical dynamics. In this application the Hartenberg- 
Denavit 4x4 transformation matrix is used to completely 
describe the positions and velocities of each point in the 


system. Appendix A develops these kinematic relationships in 
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a brief presentation, a complete review of the kinematics of 


the ERLS is available (Ref 9:pp.17-50]. 


B. KINETIC RELATIONSHIPS OF THE EQUIVALENT RIGID LINK 
SYSTEM 


The formulation of the Lagrange equations to derive the 
equations of motions is straightforward and systematic. This 
familiar approach is helpful in the complex analysis of this 
type. 

The selection of generalized coordinate systems is 
accomplished indirectly by the separation of the global 
motion into large and small motions. Formalizing this choice 
assigns the large motion gchocmmiced coordinate as theta (8), 
which represents a rotation and the small motion generalized 
coordinates (U) a column vector composed of the translations 
and rotations at the boom tip. 

The translations from the node position on the ERLS to 
the deformed position on the actual flexible member (x,y,Z) 
will become the first three small motion generalized 
coordinates. The slopes at the node or the rotational angles 
of flexible member about the axis perpendicular to the planes 
of displacement of the flexible member (U') are included as 
the last three small motion generalized coordinates. 

The transformation from the ERLS to the deformed boom 
position is accomplished by the mode ‘shape matrix. The 


derivation of the mode shape matrix is briefly presented in 
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Appendix B. The choice of nodal position and mode function 
determines the accuracy of the simulation response. 

In this representation a cubic polynomial function was 
chosen as the shape function for bending and a linear 
function for the extension and torsion. This choice provides 
the mathematical representation of our first three small 
motion generalized coordinates, and simplifies the required 
mathematics to obtain the functional representation of the 
remaining three small motion generalized coordinates. 

The system kinetic energy is the sum of the kinetic 
energies of each member of the system, this would include the 
flexible boom and the tip load (if any). The system 
potential energy results from the elastic strain energy of 
the. Elexible member, and the application of gravitational 
body ‘force (if any). Generalized forces represent the 
application of applied forces and damping forces. Appendix C 
includes the development of the equations of motion from the 
Lagrange Equations for the case of end point nodal positions 
and a polynomial shape function for six small motion 
generalized coordinates. 

The manipulations of the Lagrange equations to obtain the 
equations of motion entail extensive effort, mathematical 
dexterity, and fierce determination but the reward is a set 
of two non-linear, coupled, second-order, ordinary 


differential equations. Though both sets are coupled, one 
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describes the large motions and the other describes the small 


motions, they are presented below. 


6° +M 0 =F (2-1) 
qq qn q 
ee ae | 2 (2-2) 
nq nn n n n 
where: 
M aa = 1x 1 effective inertia matrix for large motions 
M qn = 1x 6 coupled effective inertia matrix of the 
small motion on the large motion 
F q = 1x 1 load vector for the large motion 
M na = 6x 1 coupled effective inertia matrix of the 
large motion on the small motion 
M nn = 6 x 6 effective inertia matrix for small motions 
G n = 6 X 6 gyroscopic matrix 
K n = 6 x 6 stiffness matrix 
i n = 6 x 1 load vector for the small motion 
6 = generalized coordinate of the large motion 


U = 6x 1 generalized coordinate vector for small 


motion deformations 


C. NUMERICAL INTEGRATION OF THE EQUATIONS OF MOTIONS 

The application of numerous numerical integration schemes 
is possible in the solution of the equations of motion. The 
primary consideration in the selection of the integration 
technique is the determination of the time step necessary to 


integrate the equations of motion and obtain numerical 
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stability. In this work several integrations techniques were 
tested i.e.: Runga-Kuta 5th order (both fixed and variable 
time step), Adams-Moulton (variable time step), and Gear. 
While all of these techniques achieved numerical stability 
the computational time required to obtain adequate simulation 
time was excessive. 

The formulation of the equations of motion in a manner 
similar to equations (C-~-36) and (C~-37), permits the 
application of the Sequential Integration Method [Ref. 9:pp 
79-91). In this technique an implicit integration scheme is 
employed to integrate the small motion equation, and an 
-explicit' technique is used to integrate the large motion 
equation. 

The implicit technique guarantees numerical stability in 
the case of linear equations and hence an iterative technique 
can be aves). The results of the small motion integration 
can be provided to the large motion equation which can then 
be integrated explicitly in a predictor-corrector fashion. 
The choice of time step in the case of the explicit 
integration can be governed by solution accuracy vice small 
motion frequency response. The decided advantage to this 


approach is the computational time savings. 
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III. APPLICATION OF THE ERLS TO A SPACECRAFT 
NG ELEMENT FLEXI OOM 
A. THE LOW FREQUENCY MICROWAVE RADIOMETER (LFMR) 

The low frequency microwave radiometer (LFMR) which would 
be based ona espace platform of the N-ROSS type is composed 
of a reflector, an upper reflector boom, a lower reflector 
boom, and a spacecraft boom. These booms would all be hinged 
together and owing to the requirements of space systems, 
would be developed from light weight materials. The dynamic 
modeling of such a structure has been the primary focus of 
-‘Cthe work presented here, other applications of the ERLS in 
Space technology do however exist. 

The orientations of the microwave receiving dish on board 
requires precise pointing accuracy which must be well 
~ controlled for the transmission/reception of signals. The 
accuracy of the pointing mechanism is subject to the position 
error determinations which in turn are vulnerable to the 
accelerations of orbital adjustment, antenna spinning and the 
inherent light weight materials of construction. For these 
reasons the control of vibrations is essential. 

The vibrational response or deformation resulting from 
the application of these forces must be observable to the 
controller. The choice of the small motion coordinates 
eee outlined would make these deformations available 


during start up, regular operation, and orbital station 
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keeping adjustments. Since these values are directly 
measurable this choice could simplify the controller observer 
requirements while ensuring the mission specifications are 


met. 


B. SIMULATION OBJECTIVES 

In this work the multi-bodied boom has been simplified to 
a Single member. The deformations resulting from the 
structural flexibility, have been assumed small and hence 
small angle assumptions have been implemented in the model. 
The application of body forces and internal forces have been. 
limited to gravity, and a driving torque. No damping effects 
have been implemented, nor contributions to deformations from 
thermoelastic effects upon entering or leaving the eclipse. 
The incorporation of these concerns is possible ina similar 
model but beyond the scope of the work presented here. 

The exact material of the boom would likely be a graphite 
epoxy composite with specific material properties, however 
the Bevieeion code utilizes the material considerations of 
pure aluminum to permit the intuitive interpretations 
associated with a more familiar structure material. Specific 
material considerations can be easily changed. 

The ultimate goal of this work was to complete the first 
stages of the simulation modeling. To develope the solution 


techniques and the basic simulation code (Appendices A-D) 
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which can later be expanded to eventually include multi- 


elements in the satellite model. 


C. DESCRIPTION OF THE COMPUTER SIMULATION CODE 

The IBM simulation language DSL (Dynamic Simulation 
Language) was originally chosen to model the equations of 
motions. The variety of integration schemes and FORTRAN 
compatibility were the over riding reasons’ for this choice. 
Gannon's' work in two dimensions indicated a preferred 
integration technique for the’ Adams-Moulton method (Ref. 
12@:p. 25], additionally DSL proved capable of simulating 
Gannon’s system model. 

While integration methods were not specifically the focus 
of this work, as indicated previously, the choice of 
integration technique became dependant upon the time step 
requirements for numerical stability. 

The simulation code was developed with modular FORTRAN 
subroutines then attached to a short initialization program. 
iron the three dimensional development in Appendix C the 
Lagrange equations can be formulated as equations (2-1) and 
(2-2) and then expanded in the form of equations (C-36) and 
(C-37). Where the formulation in this manner lends itself to 
the Sequential Integration Method. 


The resulting column vector: 


[6 u uu] (J=m= 
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represents the large motion acceleration 6 and the rate of 
local deformation in position and slope. Where the x y and z 
subscripts represent the local coordinates and the double x y 
and z subscripts represent the rotations about these axis. 

Each subroutine of the simulation code was developed to 
specifically determine a portion of the single element three 
dimensional analysis and due to the nature of modular 
programing the expansion of this code to a multi-body problem 
can be accommodated. The logic flow of “the coding 
(Figure 15) is provided to assist in the future development 
of multi-body simulations. 

The simulation code can be viewed from three levels; 
LEVEL I (an overview) separates the code into the primary 
focal points of INITIALIZATION, FORMULATION, INTEGRATION, 
REFORMULATION, and, OUTPUT. The reformulation is a technique 
to incorporate the inertia coupling between the large and 
small motion equations. 

LEVEL II can be viewed as the FORTRAN formulation of 
coefficients matrices, and sequential integration to solve 
for 6, 6, 6, U, U, and om 

LEVEL III is alisting of the FORTRAN utilities sub- 
programs required for the matrix formulation and manipulation 
in LEVEL II. Each of the units in LEVEL II may access the 


routines in LEVEL III. 
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D. SIMULATION RESULTS 

The computer simulation code was run with a variety of 
initial conditions these conditions were chosen to present 
the aspects of the boom flexibility and the impact this 
flexibility would have on the requirements for end point 
control. The simulation output has been presented in local 
displacement values, where these relative values represent 
the conditions at the boom tip, the point mass location. 

The simulation results presented in Appendix E are 
considered typical of results available from the simulation 
code. They represent the single large motion generalized 
coordinate and each of the six small motion generalized 
coordinates. The intent of these figures is to provide the a 
physical feeling for the plant behavior. 

The results, Figures 3-7, reflect the slowing effects on 
the large motion generalized coordinates from the point mass. 
The small motion effects, Figures 8-14, and the coupling 
between the large and small motion coordinate systems, is 
also clearly observable. Additionally the results indicate 
the effects of boom elevation angle and corresponding 


responses at four elevations, zero through 609. 
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IV. SUMMARY 


A. CONCLUSIONS 

This investigation is motivated by the needs of the 
United States Navy to be able to accurately predict the 
Meterological phenomenon within an operations theater. To 
this end the techniques currently employed for the dynamic 
analysis of large flexible space structures must be improved 
and expanded. We have presented the development and computer 
Simulation of an analysis for a flexible spacecraft boom, 
rotating in three dimensions with flexible modes in bending, 
extension and torsion. 

The nonlinear dynamic responses of the structure member 
are predicted by the numerical integration of the equations 
of motion. In this work, the spatial finite element 
discretization of the boom structure and the application of 
an assumed third order modal response provide an approximate 
solution to the equations of motion, one which lends itself 
to the problems associated with controller requirements. 

The modeling process has been applied through the 
Equivalent Rigid Link System, utilizing Lagrangian Dynamics 
to obtain the system equations of motion for numerical 
integration. The effectiveness of the integration schemes 
available for numerical integration of this type of problem 


have been briefly presented. 
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The justification for the ultimate selection ofa 
Sequential Integration Method was discussed. Areas of the 
computer code which must be developed in order to extend the 
model to multiple bodies and eventually to the controller 


design and application were also presented. 


B. RECOMMENDATIONS 
Areas which remain to be investigated utilizing the 
Simulation code include: 
1. The significance of torsion and extension in the 
total system. 
2. The signif tesnce of the way in which torque is 
applied to the boom, for example; impulse, ramp or step 


wise changes. 


3. The significance of the material properties of the 
model. . 
4. The significance of the rotational velocity. 


The single largest concept remaining to be addressed by 
this boom model is the additional members to be joined to the 
driven boom. The addition of these members’ requires 
additional application of the ERLS modeling scheme. 
Implementation requires compatibility for both members at 
the connecting joints. 

The long-term goal of this research is the design and 
control of a large flexible antenna-boom assembly for a space 


based meterological data acquisition system. Simulation 
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studies need to continue and their scope should be expanded 
to multi-body analysis. 

The control system requirements must be investigated and 
incorporated into the antenna system model. Finally a 
physical model must be developed to validate the simulation 


results. 
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APPENDIX A 


DERIVATION OF THE KINEMATICS 
REQUIRED IN THE EQUATIONS OF MOTION 


In Appendix C we presented the derivation of the 
equations of motion for the single-element, flexible boom. 
In the derivation of these equations we represent the 
position of any point in the deformed member globally by 
applying either a coordinate transformation matrix W or a 
deformation transformation matrix D and in the determination 
of the small motions the application of both transformations 
lis required. 

In order to more accurately represent the coordinate 
transformation matrix W, we will present the Hartenberg- 
Denavit transformation matrix A, and demonstrate the 


technique for obtaining the matrix W and it's derivatives. 


1 0 0 0 
a COS 8 COS 8 -SIN@ COS« SIN@ SIN« 
= a SIN 6 SIN @ COSse COSe -COSée@ SIN« 
d 0 SIN « COS « 


(A-1L) 
The theory and derivation of the Hartenberg-Denavit 
matrix formulations are defined in several of the references 


cited in Chapter One this work [Refs. 9,11-13]. The intent 
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here is to present the application of this technique. 
Variables variables used in the H-D transformation matrix are 


called joint variables and these are defined below. 


6 = the rotation of the link about the local z axis 

a = the length of the flexible link 

d= the z axis displacement between consecutive 
coordinate systems 

*® = the angular rotation between consecutive local z 


axis 


In the application of the H-D matrix, the rotation angle 
68, of the member, is always about Ene local z axis. To 
obtain this orientation we have defined two sets of local 
coordinates which will require two coordinate 
transformations. The first -transformation is a rotational 
one, we transform -:the global z axis to the eau y axis, 
(this is a rotation of « radian about the x axis). The 
transformation matrix Al represented in equation Ch 2a 


provides this rotational transformation. 


1 0 0 0 
re 0 COS 86 ~SIN6 COS«a SINS SIN« 
i 0 SIN 8 Cose COSe« -COSS SINe 
0 0 SINe COS8 


(A-2) 


eh) 


The second local coordinate system is defined at the end 
of the link, and therefore we have included the effects of 
translation. We have represented this as A , it is defined 


by equation (A-3). 


i 0 0 0 
ne a COS 8 COS 9 -~SINS COS« SINS SINe 
2 a SIN 8 SIN 9 COSS COS« -COSS8 SINe@ 
0 0 SIN « COS 9 


(A-3) 


Having defined the two H-D matrices required in the 
Simulation, we can represent the coordinate transformation 


matrix W. 


W=A,A (A-4) 


The definition of W lends itself toward the application 
of the chain rule in the determination of the derivatives of 
W. By the choice of the H-D transformation matrix we can 
easily obtain the derivatives of each A matrix as follows: 


_d 


dq A=QA48A (A-5) 


Where the derivative is taken with respect to the 
generalized coordinate, in this application q would be the 


angle theta and Q is defined by equation (A-6), 


34 


(A-6) 


Oo 

It 
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J 
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If we applied the unique conditions of the spacecraft 


rotating boom we can set the change in qd, as zero and 


eliminate the concern for order in multiplication of the 


three matrices which compose the derivative of the W matrix. 





ow 7 hs 
50, QA,A, = W, (A-7) 
and: 
aW _ 9 (A-8) 
30, 


For the time derivatives we obtain (with 89 coystant): 


, dw 9W 70 Wo 

W = —_ => — ——— «= ——.. «8 (A-9) 
ole 99, a eo, i 

meee «6G = SO n.6_A QA.@.6.A. + QA.6. A (A-10) 
dt at 31°91 "%2 te 122 a aay, 


In expression (A-10) we can consolidate the product of 


each group of terms on the extreme right hand side and 


rename them the residual acceleration ae 6? and the large 
motion acceleration Woe. In following derivations 


Subscripts for @ have been omitted since there is only one 
large motion generalized coordinate capable of motion, the 


derivatives of the other would be zero. 
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Defining the deformation transformation matrix as 


equation CAS) will accommodate three dimensional 
deformation. 
1 0 0 0 
UK] : “UK 26 Uk=5 
D = Uye> are 7 “Uy = (A=ie) 
Ue=3 Kas UK=4 i 


The partial derivative of D with respect to each of the 
small motion generalized coordinates results in the same 
matrix as the partial of the time derivative of D with 
respect to the time derivative of the small motion 
generalized coordinates, equation (A-12). In the following 
development the subscript K refers to the particular small 


motion generalized coordinate of interest. 


D = =a = 9D 


: OU, ab, 


2D 7 
7 (A-12) 
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APPENDIX B 


DERIVATION OF THE MODE SHAPE RESPONSE MATRIX @ 


The derivation of the mode shape matrix has been 
presented previously (Refs. 0971), 1 2.1. The derivation 
presented here is abbreviated and included only as a brief 
Gescription of€f the techniques involved, for complete 
derivation the reader is encouraged to review the previous 
works. 

Representing the relative position changes due _ to 
Geformations by qd we have assumed small deformations, no 
warping after torsion and all deformations u, v, w, and 8 are 
only functions of x (along the boom), the yz pane: tesmetinte 


plane after torsion and deformations only vary along the x 


direction. 
0 
u = ¥ 5 ea, ey 
d = ” i G (B-1) 
x 
Ww + Yy <i 


Applying simple beam theory, we have assumed that the 
rotary effect and the lateral shear Geflections are 


negligible. 
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(B=Ze 
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We have chosen a polynomial form to represent the shape 
function for the displacement of the single element flexible 


boom as shown below. 


Uy = ag + a,x (B-3) 

Tuy) se yee (B-4) 
2 0 1 2 3 

u. = c, + C.x +c x? + Cc x? (B-5) 
3 0 1 2 3 . 

U4 = ay (B-6) 
_ 2 Ee 
Up = b, + 2b,x + 3b4x (B-7) 
= 2 a 
Ue = Cy - 2c,x + 30.4x (B-8) 


We can now apply the boundary conditions of zero 
displacement and slope at the base where x is equal to minus 
the boom length, L, and the sign ue consistent with the 
convention for the local coordinate system. 

Substituting the boundary conditions into the _ shape 


functions gives: 
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u, (-L) = a + a,x = 0 (B-9) 
MeL) = b. + b.x + b.k° + b.x° = 0 (B-10) 
2 0 1 2 3 
u(-L) =c, + c,.x +c x2 + Cc x? = 0 (B-11) 
3 0 1 2 3 
u,(-L) = a, (B-12) 
2 
us (-L) = by + 2 box + 3 b,x (B-13) 
2 
up (Lb) = Cc) + 2 Cox + 3 Cx (B-14) 
u, (0) = a (B-15) 
un(0) = by : (B-16) 
u,(0) = Cy (B-17) 
u,(0) = ay WER Le) 
us (0) = b, : 2 (B-19) 
ue (0) = Cy (B-20) 


The coefficients a,, as, b b can now be 


2 De OGD ae 


solved for by substitution and arranged in the shape matrix 


Cc 


representation as; 


0 0 0 0 0 0 
a 0 0 0 0 0 
21 
¢ = (B-21) 
0 az5 0 az4 0 Az¢ 
O68 843 B44 345 9 
where: R = x, the variable of integration, L the length 


over which the integration will take place (the length of the 


flexible member). 
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a = 1 +R 


21 
aap OR: oR 
a = 2k Ree 
ag, = = 2 Gee) 
agg =~ ¥ (1 + R) 
ag, = LR® + 2LR* + LR 
a4, = LR° + 2LR* + LR 


For the potential energy it will be necessary to develop 


rt t t t 
u Uz, Uz-U;g, Ee accomplished by 


t 
as Pi 
differentiating the shape function equations (B-9) - (B-20) 


the expression for u 


and then solving the simultaneous equations a second time to 


produce the modified shape function matrix which is presented 


below: 
Ci 0 0 0 0 0 
: 0 0 C53 0 Cos 0 
r = 0 = 0 0 0 = (B-22) 
a2 36 
0 0 0 Cag 0 0 
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APPENDIX C 


DERIVATION OF THE EQUATIONS OF MOTION FOR THE 
SINGLE-ELEMENT, FLEXIBLE BOOM 


In derivations of this type additional definitions are 
made to Simplify the terms involved in the Lagrange 
equations. A table of terms, for both variable names and 
integral definitions (Tables I, and II respectively) have 
been included for this purpose. 

In the case of a rotating Flexible boom the large and 
small motions can be separated by application of the ERLS. 
The large motion generalized coordinate is defined as 8 and 
the aioe motion generalized coordinates as U. In this 
application U represents a column vector of six entries. We 
can imagine this vector of generalized coordinates as 
composed of two sets. The first set.ul, u2, and u3 represents 
the translational displacement of the end point in the local 
x y and z directions when deformed. The second set u4, ud, 
and u6 represents the rotation or slope of the flexed member 
at the deformed end. In accordance with this description the 
vector of generalized coordinates can be defined by equation 


(G=1)) 


6 | = ° U, YU, YU, Uys Ue uu, (C-1) 
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These generalized coordinates will be used in the 
application of the Lagrange equations to develope the 
equations of motion for both the large and small coordinate 
systems. The Lagrange equations are defined in equations 


(C-2) and (C-3). 


d [_OKE aKE aPE  _ _ 
d [aKE] _ aKE PE _ : 
dt en mo teu" —— 


Small motion generalized forces are zero in this example, 
this is consistent for a system without damping. Utilizing 
the ERLS we can define the global position of the, end point 
in terms of the local position vector (r) anda deformation 
vector (d) equation (C-4). Similarly the global position of 
the tip load is obtained by transforming a local position 
vector for the load, vector (rb), first by a deformation 
transformation matrix D and then by a coordinate 


transformation matrix W. 


R= W (nr + dq) (C-4) 
BR, = WOr, (C-5) 
In order to completely represent the coordinate 


transformation matrix we have presented the Hartenberg- 
Denavit transformation matrix A, and demonstrate the 
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technique for obtaining the matrix W and it's derivatives in 
Appendix A. Thus we continue with the development of the 
equations of motions, not encumbering the derivation of these 
equations with the derivations of others. Also presented in 
Appendix A is the definition of the deformation 
transformation matrix D. 

Continuing; we present the derivative with respect to 
time for the global positions in order to obtain the kinetic 


energy expressions. 


R= Wr +d) + Wa (Caen 
R, = WDr, + wOr, (C=7} 


e 


The kinetic and potential energy expressions equations 
(C-8) and (C-9) follow. In this development we have 
separated the large and small motion energy contributions and 


Will present them separately. 


For the large motion: 


KE=5 [fuk  Rav+ fu BR av] (C-8) 
link load 
volume volume 
PE = a fu7r*c rU dx - | u BP g av - | My. Rg av, (C-9) 
Tink link load 
length volume volume 


4 4 


Now redefining din terms of the shape matrix @ (again 
this derivation is presented in Appendix B to preserve the 


continuity here). We can rewrite equations (C-9) as: 


gh T 


= fut (ecrou) WW cctau) + cab Twiweb] + (ereu) TW web] av 


+ fulgfeorp” . "a (Dr) + (br) * W *w(Br) Js (Dr) " rs wbx) ] av, 


cC- 3 


Continuing in the ‘development we first express’ the 
derivative of the kinetic energy with respect to the time 
rate of change of the large motion coordinate, 6, equation 
(C-11), and then obtain the time rate of change of this 


expression, equation (C-12). 


QKE fuLcc+ouy *w3acc+au + (+00) Twhwad] dv + 





a6 : 
T 
Ju, [ror i 5e (Dr) $ (Dep W 54 (br, ley (C-11) 
dp 2KEy Tals T alg .yras 
at 38 = fulcad WoW(rt+OU) + (et+eU) (WoW+W W (rte) + 
(e+) "(woh (ob) + (oT (wow ced) + 
(r+) Ho wewom (60) + (rtou) (wow) (40")] av + 


i 


fu, Cebry"woacory + (ory ™(WSW + woe") (Dr) + 
(Dr) Tag 4 (br) + (Br)" Wi 5w( Br) + (Dr)" ws TyB" cr) + 


(Dr)" (We Ty + WoW) (Br) dv, (Ga) 
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Finally we express the derivative of the kinetic energy 


with respect to the large motion position as 86. 


Se =(ulsliereu asd + AM ceseuy + cody 7 ws 


o au. tae Ty 9) (a0)] 


+ (e+e) Howe wed] av + (ube [coe aga + WW, (orp + 
Hl T a eT eT e 
(Bry twow + WW) (Ded) + (Dry tao u + WW) (Ded) dv, 
(C-13) 


ry Tee ie 


OKE T A )(r+OU) + 2(r+OU) WoW (a0). + 


atl 3a] - 36 = ful (r+oU) (W 


(e+ou)Twowcad")] av + fu, [coe woe" (or + 


2(Dr) *WeW(Or) + (Dey Wow(B"r)] av, (C-14) 


The small motion equation is handled in a similar 
fashion, however it is complicated slightly by the 
deformation transformation matrix D. The product of the 
matrix terms involving D will result in scaler values, 
therefore the partial derivatives of D with respect to each 
small motion generalized coordinates will contribute 
separately. This requires the introduction of the operator ¢ 
which represents a column vector composed of the matrix 
products, where the derivatives have been taken with respect 
to each small motion generalized coordinate. 

The symbol DK represents’7 the deformation transformation 


matrix for which the partial derivative of matrix D have been 
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taken with respect to 


coordinate 


represented by K. 


equation (C-15) where K is l. 


each small motion generalized 


This is demonstrated in 


Applying the operator ¢ we can 


express the scaler terms for the derivatives of D for each K, 


1-6 by equation (C-16). 


0 0 
1 0 0 
0 0 
1 a 
D*W WD, 
pi wiwo 
TT 6 K= 
Ge, D WoW Dy = : 
| A by 
D “WWD, _ 


Separating the 
we present: 


For the link; 


3 oie 


oOo Oo AO 2 


contributions from 


Ce—15) 


o Oo Oo 2 


(E=26) 





the link and the load 


= fufe*w'wed + oTwi(r+eu)] av Tesi 
a, = fufe*i?w + wired + oTwTwod" + oTwThed + 

e7(a + wh") (r+eU)] av (C-18) 

Tt = fuferitice + ou + owed] av (C-19) 
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ae Lor - -T = fufe'w'wed’ + o*w'ted + otw'h"(c+au)] av 
U 
(C=2oR 
For the load: 
aKE, | 
aint? 
=i = Ju, [ey pew Wor, + r:D'd WDyr, J] dv, (C-21) 
K 
Tada T TT. Lams 
oe = fu fey Dy (WoW + W W)Dr, + r, Dw wD r+ 
Toten aTs TeTsT 
pr Do(W OW + WWD ee + pr D WWD re] dv, 
(C=2ay 
_ Tats TT ae z 
= fu, [epD, 8 Wor, + efor a wor.) av, (C=2en 
9KE gKE 
d L i TA ,T Tae Tahal ee TaTsT 
ael5y] ey ia fu, [e;D,W WOr, + r,D,pWWD r+ r DW WDyr, + 
K K 
TW Tact 
r, DW “WDyr, | dv, (C-24) 


Completing our development, the expressions for the 
potential energy are presented organized as the previous 


material and followed by the Lagrange equations. 


TT TTT 
5 fu uTrtcru ax + fuce + 6U) wre av fuprsD wa av, 
(C-25) 
gPE f : f TTT d 
zo u(r + $0) "Wed dv Mp D Wo dv, (C=Zep 
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Recalling equations (C-2), (C-23), and (C-26) we can make 
the appropriate substitutions and present the Lagrange 


equation for the large motion generalized coordinates. 


Fo = fu[(erourwpW"crseur + 2creeur wowed) + 


wat") dv + fu, [corp *woW" (or) + 


2(De) "wo@(Be) + (De) WowiB"r)] av, + 


(+60) "We 


el d Wa 1 8 


= fuce + $U) Woa dv = MpE,D Wea dv CC=27) 


L 


We can simplify this expression further by the 


substitution of the integral terms defined in TABLE II and 


the separation of the second time derivative of «the 
coordinate transformation matrix W by FES? force 
contributions, equation (C-28). The contributions to 


coriolis and centripetal forces are represented by Wr, termed 
the residual accelerations. The contribution to the general 


forces is represented by We. 


Ww =W + W (C-28) 


> 9 
. T T T T T T T 
F = 1,,(WgW.+ WoW) + 1,,(WoW + WOW )U + UT, (WW + WW) 
T T T Ts T Ts Tae 
+U 1,5 (WW + WW U + 22, 5(0,W)U + 2U 15 ( WW) OU + Ti 5(WWu 
T Tae TT TT T Tas 
t+ U'I,,(WyW)U + He(D WWD) + He(D WoW.D) + 2H, (D WWD) 
14 ie Tae cea TT 7 
eet. D wow ) - H,Wog - U'H,Woqg - H,D Weg (C-29) 
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Returning to the development of the small motion Lagrange 
equation the potential energy contribution remains’ to be 
presented. Again as in the previous case the contributions 
of the point load and the link will be developed separately. 
For the link: 


oe, = [ r'cru ax - fus? wig av (C-30) 


——— + = 0 (C-31) 
T aye : 


Jk ae Ts se oe T 


0 = Ip,(wwd" + ar,.whd + 1,.(wew) + TQ CWwL) + 
+ I,.(W'w_)U + I.,(W'W,)U + K,U - H.W'¢ (C-32) 
22‘ 4, 29 ac 1 2 
For the load: 
OPE 
ie TT. 7 
9KE . @PE 
a(— - TF = 0 (C-34) 
vu, ayy 
0 = H-(D-W-WD) + H, (Dy Tyo") + H, (D Tota D.) + 
5 OL" K 
T. Tes [7 
H.(D W'WO,) - H,D,W'g (C-35) 


Collecting terms and arranging the coefficients the two 
Lagrange equations can be written as the equations of motion 
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for the large and small generalized coordinates. [t is these 
equations which must be solved by the computer simulation 
code. 


For the large motion: 


T T 
1,,(WyWy) + To(W,W) + 
T | T T 
1,5(W,W,)U + U'I,,(WoW) + 
T T es T T a 7 
U T4464 G45)? + 3) + ¢ H, (D WWD)? 5 = 
T T 
et, hw )U + 
TT 
H, (D°W,W,D) 
T T T T T T 
fee Wey (Wee 0 - Ut, (WoW) - UI,,(WoW)U - 
Ta 8 T Ts s TT T Tes 
21, (W400 2U°I,,(WoW)0 - H.(D WoW D) - 2H,(D wow) + 
T T TT 
H, Wo + U HW + H,D Wo 
(e236) 
For the small motion: 
T T 
1,,(W ow) + I,,(w Ww) + 
T ee T 7 e 
ieee WU +] 8 ++ |g HL (DW wD,) 0 n 
T 1 
¢ H.(D,W WD) 
Te T 
2 1,,(W W) + 5 1,,(W W.) + Ky + ' L 
Twi ap.) ’ H.(D-wiw p F 
2¢ He(D,W WD, G He (D WW DY) 
: T T T 7 
[ Ti, (a Ww) - HAW gtg¢ H,D,W q | (C=37) 
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APPENDIX D 
LISTING OF COMPUTER SIMULATION CODE 


RRR EKER KRAARKKAARKAKAKRAARKRAAKRKEAKRAKRARAEAKAAKAAKAAARARRAAEKNREAARAEKRKRNEKRKRKKRAKAARH 


= THIS IS THE CONTROLLING PROGRAM FOR BOOM SET SDEr NES sia % 
* VARTABLES, DETERMINES THE COEFFICIENTS OK (HE fxverve.. . 
x INTEGRATION, AND CALLS THE TWO SEPARATE INTEGRATOR OF THE ai 
is LARGE AND SMALL MOTION EQUATIONS. 
. DEFINITION OF VARIABLES i! 
« x 
: BU 2) cetera e cs nee H-D JOINT VARIABLE i 
. ALPHA C2 iar H-D JOINT VARIABLE = 
AREA <3 oye sine CROSS SECTIONAL AREA OF BOOM ¥ 
ss area cy rere YOUNG'S MODULES * 
= BP SEONi. o -sae PARAMETER FOR INTEGRATION COEFFICIENTS gs 
: Pee N Meee 2 «.sPINISH TIME OF INTEGRATIOH * 
Gio oc menetene cree SHEAR MODULES . 
E GRAV Gs Sete oes GRAVITATIONAL FORCE (GLOBAL COORDINATES) i 
ai Es senetle 100 ers genre eetweme TIME STEP USED IN THE INTEGRATION is 
e 2G 8 i if: ane reer) PARAMETER FOR INTEGRATION COEFFICIENTS = 
‘ EY haere eee MOMENT OF INERTIA ABOUT THE Y AXIS * 
‘ BATA an a em ie MOMENT OF INERTIA ABOUT THE Z AXIS . 
: sy AO) ae Oaeee Rarer yc MASS DENSITY OF BOOM (KG/H**3) 7 
= beh 25 ))9Gin Geren co POINT MASS AT BOOM TIP (KG) . 
: ORES Ro Dune .... LARGE MOTION GENERALIZED COORDINATE 8 4 
‘ 0] 0) @ fA ee rare aserst » LARGE MOTION ROTATIONAL VELOCITY i 
: Lee em ent osioce/ 10sec creas INDEXING VARIABLE TO PERMIT TIMING OF PRINT * 
= 1 Gy) eee re SEC - 
e 0)74 016 ye rae nemnese LARGE MOTION ROTATIONAL ACCELERATION ~ 
. POUR hietctg ot ere tet eam LARGE MOTION FORCING FUNCTION is 
x 0 2) Ae eee SMALL MOTION AXIAL DEFLECTION * 
: C2 ee sc tee feelers SMALL MOTION TRANSVERSE DEFLECTION (Y) 5 
is CS) eine etaner eros SMALL MOTION TRANSVERSE DEFLECTION (2) i! 
. OCA e ec ee SMALL MOTION TORSIONAL DEFLECTION ie 
i LUI (2) erence rac tane car ee SMALL MOTION ROTATION ABOUT LOCAL Y AXIS * 
. GG Justa ave sree SMALL MOTION ROTATION ABOUT LOCAL Z AXIS 2 
*« x 
KKK RRR RARER EKRKEEEKEAKEAEERKEEEARKKEEKKH 


INTEGER PT,T 
REAL*8 A(2),D(2), THETA(2),ALPHA(2),AREA,E,G, 
TYY 7 24S, M0027 GRAY 
UCG A URCT U6), UZzeCTUG)7 
- Q,QDOT,Q2D0T, 
EPSLON, LOTA,H, TIME, FINTIM, TORK 


+++ + 


COMMON 
COMMON 


/BLOCK1/ A,D, THETA, ALPHA, AREA,E,G,I1YY,1I2Z, MU, GRAV 
/BLK15/ Q,QDOT,Q2DO0T 
COMMON /BLK16/ U,UDOT,U2DOT 
GRAV=0.DO0 
TORK=100.D0 
: MU(2)=44.1 
MU(2)=0.D0 
PINTIH =2. 00 
H = .0006D0 
PT=34 
EPSLON 


25D0 
IOTA = 0 


= 0. 
0.5D 
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10 


A(1)=0.0 

D(1)=0. 

THETA(1)=Q 
ALPHA(1)=2. *ATAN(1. ) 


A(2)=5.486 
D(2)=0. 
ALPHA(2)=0. 


AREA=1.9E-3 
E=70.E9 
G=26.0E9 
IYY=34.63E-6 
IZZ=34.63E-6 
MU(1)=2710.D0 


Q = 0.D0 

QDOT = 0.D0 

Q2D0T = 0.D0 

DO 5 I=1,6 
U(I)=0.D0 
UDOT(I)=0.D0 
U2DOT(I)=0.D0 


THETA(2)=0.D0 
T=0 
DO 10 TIME=0,FINTIM,H 
T=T+l1 
CALL CCOEF(EPSLON,IOTA,H) 
CALL INTSML( TIME, TORK) 
CALL INTLRG(TIME,H, TORK) ’ 
PeGllLHE. EO 40.00 ent RINTL( TIME) 
IF(T.EQ. PT) THEN 
CALL RINT1( TIME) 
T=0 
ENDIF 
CONTINUE 


END 
SUBROUTINE INTSML( TIME, TORK) 


Pe ee a a a! 


* 
* 
* 
* 
* 
* 
* 
x 
x 
* 
* 
® 
*« 
x 
* 
* 
* 
* 
* 
* 
* 
* 


THIS SUBROUTINE IMPLICITLY INTEGRATES THE SMALL MOTION EQUATION 


ADDITIONAL DEFINITION OF VARIABLES 


Pees ss seapeene. 6 THE COEFFICIENTS FOR NUMERICAL INTEGRATION 

Fleece meer ae eu. THE SMALL MOTION FORCE VECTOR 

PON Nieeeetecs apes ome REFORMED SMALL MOTION FORCE VECTOR 

|HOlovm Gi6he Oho OA ae THE LARGE MOTION FORCE VECTOR 

NeeMears mote ea ciaiss. oie « GYROSCOPIC MATRIX 

(1 ea oe STIFFNESS MATRIX 

Lo I Ne ene coe REFORMED SMALL MOTION INERTIA AND COUPLING 
MATRIX 

PUNINicwemeworsc: sinter en, SMALL MOTION INERTIA MATRIX 

pl Ol ons Gunns 5 5H COUPLING INERTA MATRIX 

FU Ghar oa 6. et autacaire LARGE MOTION INERTIA MATRIX 

Oe Petes) eile se) 6 seine TEMPORATY MATRICES FOR INTERMEDIATE VALUES 

SNS D otoy Gl Saree eeey ree DUMMY ARGUMENT USED IN TROUBLE SHOOTING 

WKAREA....55.%. STORAGE FOR IMSL LEQT2F ROUTINE 


(USED FOR LINEAR EQUATION SOLVING) 


wR TRH R AHR UEKKHKK HK KH 
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* %* *© + + © © # + 4 © 4 © & & & HB HH H 


10 


20 


30 


40 


REAL*8 MQQ,MQN(6),MNN(6,6),GN(6,6),KN(6,6),FN(6),FQ, 


MN(6,6),FNN(6), 
A0,A1,A2,A3,A4, 


+++ ++ 


COMMON /BLK16/ U,UDOT, 


U(6),UDO0T(6),U2D0T(6), 


A5,A6,A7, 


U2DeT 


WKAREA(54),PRD(6,6),TIME, TORK, 
TEMP1(6),TEMP2(6), TEMP3(6), TEMP4(6), TEMP5(6. 56) 


COMMON /BLOCKS/ MQQ,MQN,MNN,GN,KN,FN,FQ 


COMMON /BLOCK7/ MN, FNN 


COMMON /BLOCK8/ AO,A1,A2,A3,A4,A5,A6,A7 


CALL TRANSF(TIME) 
CALL WMATRX 


CALL WDERIV 


IF (TIME.EQ.0.D0)CALL LRGCOF( TORK ) 


CALL SMLCOF 
CALL REFORM. 
DO 10 f#1,6 


TEMP1(T) 
TEMP ZGr) 


BL UCL) 


AO*U(L) + A2*UDOT(L) + A3*U2DOT(I) 
+ A4*UDOT(L) + AS*U2D0T(L) 


CALL MATMUL(PRD,MN,6,6,TEMP1,6,1, TEMP 3) 


CALL MATMUL(PRD,GN,6,6, TEMP2,6,1,TEMP4) 


DO 20 I=1,6 
TEMP1(T1) 
TEME ets) 


UCT) 


DO 20 J=1,6 
TEMPS(I,J) = K 


i) 
N = 6 
IDGT = 1 


M 


FNN(I) + TEMP3(1) 


N(I,J) + AO*MN(I,J) 


+ TEMP4(1) 


CALL LEQT2F( TEMPS,.M,N.N, TEMP1,IDGT,WKAREA, IER) 


DO 30 I=1,6 
UCE) = TEMNPI(1L) 


DO 40 I=1,6 
TEMP3(T1) 
U2ZDO EG) 


UZDOTCI) 
AO*(U(T) 


UDOT(I) = UDOT(I) + A6*TEMP3(T) 


END 


— tone 2 Glo 
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- A2*UDOT(I) 


+ A7*U2D0T(TI) 


+ A1L*GN(I,J) 


- AS*TEMP3(1) 


SUBROUTINE INTLRG (TIME,H, TORK) 


RRMRERKRKKEKKEKRKEKEREKEKRKERKEKEKEKKEKAKKEKEKKKE KKK KEKE KRKKK KK 


x 
x 
x 
x 
x 
* 
* 
x 
* 
* 
* 
x 
x 
x 
” 
x 
x 
* 
* 
® 


+++ + 
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THIS SUBROUTINE INTEGRATES THE LARGE MOTION EQUATION EXPLICITLY 


ADDITIONAL DEFINITION OF VARIABLES 


[slo Panera hts Se THE COEFFICIENTS FOR NUMERICAL INTEGRATION 

EN ss cissspletto me tees .THE SMALL MOTION FORCE VECTOR 

PQS ues oe ew THE LARGE MOTION FORCE VECTOR 

O(a = REFORMED LARGE MOTION FORCE 

GNo »+.-GYROSCOPIC MATRIX 

I Nicmeestoa ss 6 15s STIFFNESS MATRIX 

eI )6 6 1G SiGncemenen 5c REFORMED SMALL MOTION INERTIA AND COUPLING 
MATRIX 

MNN..........-SMALL MOTION INERTIA MATRIX 

PAGIN os eons oie « COUPLING INERTA MATRIX 

PG Oraet. secre LARGE MOTION INERTIA MATRIX 

8) 2) Re Go oc eeo TEMPORATY MATRIX FOR INTERMEDIATE VALUES 

PRD... 3. , > DUMMY ARGUMENT USED IN TROUBLE SHOOTING 


RRR AKER KKK RRR RHR KKK RK KKRKEKAMKEKKKERAKRKEKAK K 


REAL*8 MQQ,MQN(6),MNN(6,6),GN(6,6),KN(6,6),FN(6),FQ, 


U(6),UD0T(6),U2D0T(6), 
Q,QD0T,Q2D0T, 
A0,A1,A2,A3,A4,A5,A6,A7, 
H,TEMP1,PRD(6,1),FQN, TORK, TIME 


COMMON /BLOCKS/ MQQ,MQN,MNN,GN,KN,FN,FQ 
COMMON /BLK15/ Q,QDOT,Q2D0T : 

COMMON /BLK16/ U,UDOT,U2DOT 

COMMON /BLOCK8/ AO,A1,A2,A3,A4,A5,A6,A7 


Q= Q + H*QDOT + A3*Q2D0T/AO 
QDOT = QDOT + A6&*Q2D0T 


CALL TRANSF (TIME) 
CALL WMATRX 
CALL WDERIV 
CALL LRGCOF( TORK) 


TEMP1=0.D0 
DO 10 I[=1,6 
TEMP1 = TEMP1 + MQN(I)*U2DOT(I) 


FQN = FQ - TEMP1 
Q2D0T = (1.D0/MQQ) *FQN 


QDOT = QDOT + Q2D0T*A7 
Q = Q + Q2D0T/A0 


END 


ah, 


x 
* 
* 
* 
* 
*% 
x 
x 
x 
x 
x 
* 
* 
x 
* 
* 
x 
x 
bd 
x 


SUBROUTINE LRGCOF( TORK) 


RRMA AKRRARARRARARKRARAARARACRKRRARKRKERKTRARRRARARARRAARARARKERKRRAKRNRRKRRRKRRKRRKKRKRARK 


x xX 
* THIS SUBROUTINE CALCULATES VALUES FOR COEFFICIENT MATRICES FOR * 
* LINK ONE. * 
x bi 
* ADDITIONAL DEFINITION OF VARIABLES * 
x x 
* Kl -7 . 23 THE COEFFICIENTS FOR NUMERICAL INTEGRATION * 
* CONST... ere SINGLE VALUED TERM USED IN SUMMATION FOR * 
* FQ VALUE * 
* BG ha eee DERIVATIVE OF DFORM MATRIX WRT SMALL MOTION * 
‘ COORDINATE OF INTEREST * 
* DCT........... TRANSPOSE OF DC MATRIX * 
* COLO ek se THE RATE OF DEFORMATION TRANSFORMATION MATRIX®* 
* BOOTT an. one THE TRANSPOSE OF THE DDOT MATRIX * 
‘ DFORM......... THE DEFORMATION TRANSFORMATION MATRIX * 
* DFORMT..... ...THE TRANSPOSE OF THE DFORM * 
* oh re re THE SMALL MOTION FORCE VECTOR * 
* Qn ehe ee THE LARGE MOTION FORCE VECTOR * 
* RON i arietyateacls REFORMED LARGE MOTION FORCE * 
* Neat ereetote,. GYROSCOPIC MATRIX * 
* KN aeyenens see STIFFNESS MATRIX * 
* MATS" 4. oe VARIOUS DIMENSIONED MATRICES OF INTERMEDIATE* 

* VALUES (***) FOR THE VARIOUS SIZES * 
* MEA ee REFORMED SMALL MOTION INERTIA AND COUPLING * 
: MATRIX * 
* MNIN Gises a tigciv sida SMALL MOTION INERTIA MATRIX * 
* MOND. aise nee COUPLING INERTA MATRIX * 
* MOG face ahs ..LARGE MOTION INERTIA MATRIX * 
* TEMP1....:....TEMPORATY MATRIX FOR INTERMEDIATE VALUES * 
* ONE rane crt aee ARRAY OF UNITY NECESSARY IN OBTAINING * 
* DERIVATIVE OF DFRM MATRIX = 
* BRD Shh cete ae DUMMY ARGUMENT USED IN TROUBLE SHOOTING * 
* Otani tee HARTENBERG DENAVIT TRANSFORMATION MATRIX * 
* WDM Elec ace ck or winan DERIVATIVE OF W WRT @ * 
* DOT yaoi DERIVATIVE OF W WRT TIME * 
* WRESID........RESIDUAL ACCELERATION TRANSFORMATION * 
: BY WA ict cuoca ae VARIOUS KERNELS FOR THE BOOM VOLUME INTEGRAL * 
* Mla eaten. LOWER LIMIT OF BOOM VOLUME INTEGRATION * 
' KO esau UPPER LIMIT OF BOOM VOLUME INTEGRATION * 
* oy ee DUMMY ARGUMENTS * 
x x 
x x 


RRARARRRRRRRAAAAARARAT ARATE RARER KAKA 


INTEGER INDEX,NR,NC,NR2, 


NC2,TERM 


REAL*8 A(2),D(2),THETA(2),ALPHA(2),AREA,E,G, 


UT(6), 
Q,QDOT,Q2D0T, 
W(4,4), 


++ eeterere ++ ret 


WPRM(4,4),WDOT(4, 
WD1(4,4),WD2(4,4),WD3(4,4),WD4(4,4),WD6(4,4), 
WD9(4,4),WD10(4,4),WD13(4,4),WD14(4,4), 
MQQ,MQN(6),MNN(6, 
MAT(6),MAT4A(4),MAT4B(4),MAT64A(6,4),MAT64B(6,4),MATPRD(6,6), 
CONST, 2(6,6),Y¥(6, 
PTMASS(4),GRAVEC(4),WT(4,4),WPRMT(4,4),PRD(6,6),0C(4,4), 
DFRM(4,4),DFRMT(4,4),DDOT(4,4),DDOTT(4,4),DCT(4,4), TORK 


I¥Y,12Z2,MU(2),GRAV, 
U(6),UDOT(6),U2D0T(6), 


4), WRESID(4,4), 


6),GN(6,6),KN(6,6),FN(6),FQ, 


6),XL,XU,ONE(6),MAT44A(4,4),MAT44B(4,4), 
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COMMON /BLOCK1/ A,D, THETA, ALPHA, AREA,E,G,IYY,1Z2Z,MU,GRAV 
COMMON /BLK15/ Q,QDOT,Q2D0T 

COMMON /BLK16/ U,UDOT,U2D0T 

COMMON /BLOCK3/ W 

COMMON /BLOCK4/ WD1,WD2,WD3,WD4,WD6,WD9,WD10,WD13,WD14 
COMMON /BLOCK5/ MQQ,MQN,MNN,GN,KN,FN,FQ 

COMMON /BLOCK6/ WPRM,WDOT,WRESID 


XL=-A(2) 
XU=0.D0 


MQQ=0.0D0 
FQ=TORK 
DO 3 I=1,6 

3 MQN(I)=0.D0 


DO 5 J=1,4 
GRAVEC(J)=0.D0 
5 PTMASS(J)=0.00 


GRAVEC(4)=GRAV 
PTMASS(1)=MU( 2) 
TERM=0 


DO 7 J=1,6 
7 ONE(J)=1.D0 


CALL DFORM(1,7,U,DFRM) 

CALL DFORM(0,7,UDOT,DDOT) 

CALL TRANSP( TERM,OFRM,4,4,DFRMT) 
CALL TRANSP( TERM,DDOT,4,4,DDO0TT) 


< CALL TRANSP(TERM,U,6,1,UT) 
CALL TRANSP(TERM,W,4,4,WT) 
CALL TRANSP( TERM, WPRM,4,4,WPRMT) 


RRM RRR RA RRA RRR ARRRRARRRRARARRAERARRRRRRRAARRKKR 


* CREATING THE COEFFICIENT MATRIX MQQ 


i220 20 ff 222 eee eee ee eee eee eee RR E SE RES ESE SS ESESSESESSSSLAOSSCOS eS See ee ee 


INDEX = 11 
CALL BINT( TERM, XL,XU,Z,Y,1,1,INDEX,WD6,CONST) 


MQQ=MQQ+CONST 


INDEX = 12 

CALL BINT(TERM,XL,XU,Z,Y,1,6, INDEX, WD6,MAT) 
CALL MATMUL(PRD,MAT,1,6,U,6,1,CONST) 
MQQ=MQQ+CONST 


INDEX = 21 

CALL BINT(TERM,XL,XU,Z,Y,6,1,INDEX,WD6,MAT) 
CALL MATMUL(PRD,UT,1,6,MAT,6,1,CONST) 
MQQ=MQQ+CONST 


INDEX = 22 

CALL BINT( TERM, XL,XU,2,Y.6,6, INDEX,WD6,MATPRD) 
CALL MATMUL(PRD,MATPRD,6,6,U,6,1,MAT) 

CALL MATMUL(PRD,UT,1,6,MAT,6,1,CONST) 
MQQ=MQQ+CONST 


* THE CONTRIBUTION OF THE POINT LOAD 


CALL MATMUL(PRO,WD6,4,4,DFRM,4,4,MAT44A) 
CALL MATMUL(PRD,DFRMT,4,4,MAT44A,4,4,MAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY [TS TRANSPOSE. 
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CONST=MU(2) *MAT44B(1,1) 
MQQ=MQQ+CONST 


RRRRRRRRAAAARAEAARARAERAREARRRRARARARAAKRTRTAAERRARERTRAEKRRRAKRERKAEARARERAERRARERERER 


® 


CREATING THE COEFFICIENT MATRIX HQN 


RRR ARRRAARRRKRARAREEKARRAAAEAATAAERAARRAAEEKRRKRKKRAKKEAEAERRRRRAEAEERER 


x 


9 


INDEX=12 


CALL BINT(TERM,XL,XU,2,Y,1,6,INDEX,WD2, MAT) 
DO 8 I=1,6 
MQN(I)=MQN(I) + MAT(I) 


INDEX = 22 
CALL BINT(TERM,XL,XU,2,Y,6,6, INDEX,WD2,MATPRD) 
CALL MATMUL(PRD,UT,1,6,MATPRD,6,6,MAT)- 
DO 9 I=1,6 
MON(I) = MQN(I) + MAT(T) 


THE CONTRIBUTION OF THE POINT LOAD 


DO 10 K=1,6 
CALL DFORM(0,K,ONE,DC) 
CALL MATMUL(PRD,WD2,4,4,DC,4,4,MAT44A) 
CALL MATMUL(PRD, DFRMT,4,4,HAT44A,4,4, MAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 


10 


MQN(K)=MQN(K) + MU(2)*MAT44B(1,1) 


RREACRARAKAEARAARARRACRARAAERRAARRAARARARARKRRKRAAKRRAKCKRRARRARAKAARERARKARKERAKRAEKKARRTARKAAKRK 


x 


CREATING THE COEFFICIENT MATRIX FQ 


RRERWMARRREAKTATERAKAEEKAEATRAEARAKEEEKATRERTATAARRAACACRAERKERAAEARERTACARNTRAREARRAEKAKRARAARKANK 


x 


INDEX = 11 
CALL BINT( TERM, XL, XU,2Z,Y,1,1, INDEX, #D14,CONST) 
FQ=FQ-CONST 


INDEX = 12 

CALL BINT(TERUM, XL, XU,2,Y,1,6, INDEX, WD10,MAT) 
CALL MATMUL(PRD,MAT,1,6,UD0T,6,1,CONST) 
FQ=FOQ=2. *CONST 


CALL BINT(TERM, XL,XU,2,Y,1,6, INDEX, WD14,MAT) 
CALL MATMUL(PRD,MAT,1,6,U,6,1,CONST) 
FQ=FQ-CONST 


INDEX = 21 

CALL BINT( TERM, XL,XU,Z,¥,6,1, INDEX, WD14, MAT) 
CALL MATMUL(PRD,UT,1,6,MAT,6,1,CONST) 
FQ=FQ-CONST 


INDEX = 22 

CALL BINT( TERM, XL, XU,2,Y¥,6,6,INDEX, WD10,MATPRD) 
CALL MATMUL( PRD, MATPRD,6,6,UDOT,6,1,MAT) 

CALL MATMUL(PRD,UT,1,6,MAT,6,1,CONST) 

FQ-FQ-2. CONST 


CALL BINT(TERM,XL,XU,2,Y¥,6,6, INDEX, W014, MATPRD) 
CALL MATMUL(PRD,MATPRD,6,6,U,6,1,MAT) 

CALL MATMUL(PRD,UT,1,6,MAT,6,1,CONST) 
FQ=FQ-CONST 


THE CONTRIBUTION OF THE POINT LOAD 
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CALL MATMUL(PRD,WD14,4,4,DFRM,4,4,MAT44A) 
CALL MATMUL( PRD, DFRMT,4,4,MAT44A,4,4,MAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 


CONST=MU(2) *MAT44B(1,1) 
FQ=FQ-CONST 


CALL MATMUL(PRD,WD10,4,4,DDOT,.4,4,MAT44A) 
CALL MATMUL( PRD, DFRMT, 4,4,MAT44A,4,4,MAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 


CONST=2.D0*MU( 2) *MAT44B(1,1) 
FQ=FQ-CONST 


* THE CONTRIBUTION OF POTENTIAL ENERGY 
INDEX = 10 
* Hl, WD1l WILL BE CARRIED AS A DUMMY ARGUMENT 
CALL BINT( TERM, XL,XU,2,¥,1,4, INDEX, WD1,MAT4A) 
CALL MATMUL(PRD,MAT4A, 1,4, WPRMT, 4, 4,MAT4B) 
CALL MATMUL( PRD, MAT4B,1,4,GRAVEC, 4, 1,CONST) 
FQ=FQ+CONST 
INDEX = 20 
* H2, WD1l WILL BE CARRIED AS A DUMMY ARGUMENT 
CALL BINT( TERM, XL, XU,2Z,Y,6,4, INDEX, WD1,MAT64A) 
CALL MATMUL(PRD,UT,1,6,MAT64A, 6,4, MAT4A) 
CALL MATMUL(PRD,MAT4A, 1,4, WPRMT,4,4,MAT4B) 
“CALL MATMUL(PRD,MAT4B, 1, 4,GRAVEC, 4, 1,CONST) 
FQ=FQ+CONST 
* THE CONTRIBUTION OF THE POINT LOAD TO THE POTENTIAL ENERGY 
CALL MATMUL( PRD, DFRMT, 4,4, WPRMT, 4,4, MAT44A) 
CALL MATMUL( PRD, MAT44A, 4, 4,GRAVEC, 4,1,MAT4B) 


CALL MATMUL(PRD, PTMASS,1,4,MAT4B,4,1,CONST) 
FQ=FQ+CONST 


RETURN 
END 
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SUBROUTINE SMLCOF 


Le 2 eee eee eee eee EES ES RSS SRE SSE EEE ERE RRS EERE ERR SRE SESE SESE ERS RRS S REE EES BS @ 


x x 
* THIS SUBROUTINE CALCULATES VALUES FOR COEFFICIENT MATRICES FOR * 
* LINK ONE. * 
x nx 
* ADDITIONAL DEFINITION OF VARIABLES ‘ 
x Rn 
* eae ere eye THE COEFFICIENTS FOR NUMERICAL INTEGRATION * 
* CONST.........SINGLE VALUED TERM USED IN SUMMATION FOR ‘ 
‘ FQ VALUE * 
* DC............DERIVATIVE OF DFORM MATRIX WRT SMALL MOTION * 
: COORDINATE OF INTEREST * 
* Det Aan oe TRANSPOSE OF DC MATRIX * 
* DOCie ay. ....THE RATE OF DEFORMATION TRANSFORMATION MATRIX* 
* DOOR Thee ...THE TRANSPOSE OF THE DDOT MATRIX * 
* DEORE wi-oac anaes THE DEFORMATION TRANSFORMATION MATRIX * 
* DFORMT. eon THE TRANSPOSE OF THE DFORM * 
* A org ene ee THE SMALL MOTION FORCE VECTOR * 
* FO. i0i: ......THE LARGE MOTION FORCE VECTOR * 
* BONER sass «sce REFORMED LARGE MOTION FORCE * 
* GR aries GYROSCOPIC MATRIX * 
* KN cece STIFFNESS MATRIX 7 
‘ HAT See. eee VARIOUS DIMENSIONED MATRICES OF INTERMEDIATE * 
* VALUES (***) FOR THE VARIOUS SIZES * 
* MN............REFORMED SMALL MOTION INERTIA AND COUPLING * 
* MATRIX * 
* HNN; sa eee SMALL MOTION INERTIA MATRIX * 
‘ HON. cgis soe COUPLING INERTA MATRIX * 
* HOO nts) aca ee LARGE MOTION INERTIA MATRIX * 
* TEMP Tar ee TEMPORATY MATRIX FOR INTERMEDIATE VALUES ‘ 
* ONE os eGeaeee ARRAY OF UNITY NECESSARY IN OBTAINING * 
* DERIVATIVE OF OFRM MATRIX , " 
* PROS ee DUMMY ARGUMENT USED IN TROUBLE SHOOTING * 
* Pg tce earn A. HARTENBERG DENAVIT TRANSFORMATION MATRIX ‘ 
‘ @PRM Sg. oc +, DERIVATIVE OF W WRT 0 * 
‘ WOOT Ac... ce DERIVATIVE OF W WRT TIME * 
* WREGID.< 5%: RESIDUAL ACCELERATION TRANSFORMATION * 
‘ WI=T4. . oonee VARIOUS KERNELS FOR THE BOOM VOLUME INTEGRAL * 
* {a eee: LOWER LIMIT OF BOOM VOLUME INTEGRATION * 
* eee eee UPPER LIMIT OF BOOM VOLUME INTEGRATION * 
* VAC Me hia ee DUMMY ARGUMENTS ‘ 
x w x 
*® « 


RRAAKKKKKRKKKREKRARKKKKRKAEKARRAEARRAKRKARKRKKRAEAEKKARKKRAERAKKKRKRAKRRAARAKRARRRAKARHKRRRAN 


INTEGER INDEX,NR,NC,NR2,NC2, TERM 
REAL*8 A(2),0(2),THETA(2),ALPHA(2),AREA,E,G, 
TYY,1Z2Z,MU(2),GRAV, 
w(4,4), 
U(6),UD0T(6),U2D0T(6), 
UT(6), 
Q,.QDOT,Q2D0T, 


WPRM(4,4),WDOT(4,4),WRESID( 4,4), 
WO104,4),W02(4,4),WD03(4,4),WD4(4,4),WD6(4,4), 
WO9(4,4),WD010(4,4),WD13(4,4),WD14(4,4), 
MQQ.MQN(6),MNN(6,6),GN(6,6),KN(6,6),FN(6),FQ, 
MAT(6),MAT4A(4),MAT4B(4),MAT64A(6,4),MAT64B(6,4),MATPRD(6,6), 
CONST, 2(6,6),¥(6,6),XL,XU,ONE(6),MAT44A(4,4),MAT44B( 4,4), 
PTMASS(4),GRAVEC(4),WT(4,4),WPRMT(4,4),PRD(6,6),0C(4,4), 
DFRM(4,4),DFRMT(4,4),DD0T(4,4),DD0TT(4,4),0CT(4,4) 


++ er etrtretrere +e eee t 


COMMON /BLOCK1/ A,D,THETA, ALPHA, AREA,E,G,IYY.I1ZZ,MU,GRAV 
COMMON /BLK15/ Q,QDOT,Q2D0T 

COMMON /BLK16/ U,UDOT,U2D0T 

COMMON /BLOCK3/ W 

COMMON /BLOCK4/ WD1,WD2,WD3,WD4,WD6,WD9,WD10,WD13,WD14 
COMMON /BLOCKS/ MQQ,MQN,MNN,GN,KN,FN,FQ 

COMMON /BLOCK6/ WPRM, WOOT, WRESID 
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XL=~A(2) 
XU=0.D0 
DO 3 [=1,6 
FN(I)=0.D0 
DO 3 J=1,6 
GN(I,J)=0.D0 
KN(I,J)=0.D0 


3 MNN(I,J)=0.D0 
DO 5 J=1,4 
GRAVEC(J)=0.D0 
2 PTMASS(J)=0.D0 


GRAVEC (4)=GRAV 
PTMASS(1) =MU( 2) 
TERM=0 
DO 7 J=1,6 
7 ONE(J)=1.D0 
CALL DFORM(1,7,U,D0FRM) 
CALL DFORM(0,7,UDOT,DDOT) 
CALL TRANSP( TERM,DFRM, 4,4,DFRMT) 
CALL TRANSP(TERM,DDOT,4,4,DDOTT) - 
CALL TRANSP(TERM,U,6,1,UT) 
CALL TRANSP( TERM,W,4,4,WT) 
CALL TRANSP( TERM, WPRM, 4,4, WPRMT) 


RAAAAARARARARARAAARAAKKRARAEKRAARRREEARKRRAAKARARRARRERARRARARRRERERRARR 


* CREATING THE COEFFICIENT MATRIX MNN 
RARKAAAAAAAAAAEARAAKRAAAARERAE REAR RAREREREA RARER RES 
INDEX = 22 ; 
CALL BINT( TERM, XL,XU,2Z,Y,6,6, INDEX,WDO1,MATPRD) 
DO 12 I=1,6 -, 
DO 12 J=1,6 
12 MNN(I,J)=MNN(I,J) + MATPRD(I,J) 


* THE CONTRIBUTION OF THE POINT LOAD 


DO 15 K=1,6 
CALL DFORM(0,K,ONE,DC) 
CALL TRANSP( TERM,DC,4,4,DCT) 
CALL MATMUL(PRO,DCT,4,4,W8D1,4,4,MAT44A) 
CALL DFORM(0,K,ONE,DC) 
CALL MATMUL(PRD,MAT44A,4,4,DC,4,4,HAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 


15 MNN(K,K)=MNN(K,K) + MU(2) *MAT44B(1,1) 


RARAAARARAKRAARARARARAKRAARARAAARARARARAKRARARKRAARAAKRRARKRAAKRAARRARARAKRARAKRAARAKRARAKRARARARARARARKRARN 


* CREATING THE COEFFICIENT MATRIX GN (UDOT) 


RARAAARARARAAKRARAKRAAAARAARAAARARAARARAKRARARKRARARAKRARAKRARARARARAAKRAAKRAKRAARARRARARAARARRARARR 


INDEX = 22 
CALL BINT( TERM, XL,XU,2,Y,6,6, INDEX, WD9,MATPRD ) 
DO 22 [=1,6 
DO 22 J=1,6 
22 GN(I,J)=GN(I,J)+2*MATPRD(I,J) 


* THE CONTRIBUTION FROM THE POINT LOAD 


DO 30 J=1,6 
CALL DFORM(0,J,ONE,DC) 
CALL MATMUL(PRD,WD9,4,4,DC,4,4,MAT44A) 
CALL TRANSP(TERM,DC,4,4,DCT) 
CALL MATMUL(PRD,DCT,4,4,MAT44A,4,4,MAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 
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CONST=2*MU( 2) *MAT44B(1,1) 
30 GN(J,J)=GN(J,J)+CONST 


RRERRRARERARERERRERRERERRAREREERERERERRARAEARRAEAREARRARRERAARRAEARARAAAS 


* CREATING THE COEFFICIENT MATRIX KN(U) 


RRA RARARRRRERRERARERRRRAREREREERER ARERR ERAERRERRARRAARARARARRERARAERAARRR 


INDEX = 22 
CALL BINT( TERM, XL,XU,Z,Y¥,6,6, INDEX, #D13,MATPRD) 
DO 24 [=1,6 
DO 24 J=1,6 
24 KN( I,J) =KN(I,J)+MATPRD(I,J) 


* THE CONTRIBUTION OF POTENTIAL ENERGY FROM THE LINK 
* K1, WD1 WILL BE CARRIED AS A DUMMY ARGUMENT 


INDEX = 1 
CALL BINT( TERM, XL, XU,2,Y,6,6, INDEX, WD1,MATPRD) 


DO 28 [=1,6 


DO 28 J=1,6 
28 KN(I,J)=KN(I,J) +MATPRD(I,J) 


* THE CONTRIBUTION FROM THE POINT LOAD 


DO 32 J=1,6 
CALL DFORM(0,J,ONE,DC) - 
CALL MATMUL(PRD,WD13,4,4,DC,4,4,MAT44A) 
CALL DFORM(0,J,ONE,0C) 
CALL TRANSP(TERM,DC,4,4,0CT) 
CALL MATMUL(PRD,DCT,4,4,MAT44A,4,4,MAT44B) 


* ONLY ELEMENT MAT44B(1,1) WILL SURVIVE FOLLOWING POST MULTIPLICATION BY 
* LOCAL POSITION VECTOR AND PRE MULTIPLICATION BY ITS TRANSPOSE. 


CONST=MU( 2) *MAT44B(1,1) 
32 KN(J,J)=KN(J,J)+CONST 


RRARARAARAKAAAAARARKAAARAARAARAARKARARARARAARARARARARARARARARARARAARARARAEARAARAKAREARKARNR 


* CREATING THE COEFFICIENT MATRIX FN 


RAEAAAAAARARAARARARKRARARARARARARAARARARARARARARRARARAARARAARARARARAARAKRARARRARRAAKRAARARR 


INDEX = 21 
CALL BINT( TERM, XL,XU,2Z.,Y,6,1,INDEX,WD13,MAT) 
DO 20 J=1,6 

20 FN(J)=FN(J)-MAT(J) 


* THE CONTRIBUTION OF POTENTIAL ENERGY FROM THE LINK 
* H2, WD1 WILL BE CARRIED AS A DUMMY ARGUMENT 
INDEX = 20 


CALL BINT(TERM,XL,XU,2,Y,6,4, INDEX, WD1,MAT64A) 
CALL MATMUL(PRD,MAT64A,6,4,W0T,4,4,MAT64B) 
CALL MATMUL(PRD,MAT64B,6,4,GRAVEC,4,1,MAT) 
DO 26 J=1,6 
26 FN(J)=FN(J)+MAT(J) 


“ THE CONTRIBUTION OF THE POINT LOAD TO THE POTENTIAL ENERGY 


DO 36 J=1,6 
CALL DFORM(0,J,ONE,DC) 
CALL TRANSP(TERM,DC,4,4,DCT) 
CALL MATMUL(PRD,DCT, 4,4,WT,4,4,MAT44A) 
CALL MATMUL(PRD, MAT44A,4,4,GRAVEC,4,1,MAT4A) 
CALL MATMUL(PRD, PTMASS,1,4,MAT4A,4,1,CONST) 
36 FN(J)=FN(J)+CONST 


RETURN 
END 
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SUBROUTINE WMATRX 


RARRAAKRAEEAKEKKRARARKRAKRAKRAARKRRRRARRARAARKRARARRAKARAKRRAAERRKRARRRKRARAAARAKRKRERRRER 


THIS SUBROUTINE MULTIPLIES THE HARTENBERG-DENAVIT MATRICES 
TOGETHER TO OBTAIN THE COORDINATE TRANSFORMATION MATRIX 


HD...........-+THE HARTENBERG DENAVIT TRANSFORMATION MATRIX 


& ® 
& ® 
& x 
* x 
S ADDITIONAL DEFINITION OF VARIABLES * 
R * 
& ® 
a veeeeesee TRANSFORMATION MATRIX * 
x ® 
® x 


RRRAARKKRKARRARKRKRRARKRAKRAARARAEARARARAKRRAARRAARRAAERARRARARKRARRARRARRRARRRAKRKRAKRE 


REAL*8 8(4,4),HD(2,4,4),8SUM(4,4),SUM 


COMMON /BLOCK2/ HD 
COMMON /BLOCK3/ W 


DO 10 L=1,4 


DO 10 M=1,4 
WSUM(L,M)=0.D0 
10 W(L.M)=0.D0 
DO 40 J=1,4 
DO 30 K=1,4 
SUM=0.0 
DO 20 M=1,4 
SUM = SUM + HD(1,J,M) * HD(2,M,K) 
20 CONTINUE 
WSUM(J,K)=SUM 
30 CONTINUE 
40 CONTINUE 
DO 50 I=1,4 
DO 50 J=1,4 
50 W(I,J)=WSUM(I,J) 
END 


SUBROUTINE WDERIV 


RRAAKRAKEARAKAKEKERAKERRAARRARKRKRARAAARARARAAERKRKRKEAREKEKRAKRARRARAAERAKRARAKNKAEAHEKEHK 


x : x 
a THIS SUBROUTINE CREATES THE DERIVATIVES OF THE COORDINATE MATRIX * 
. W WRT Q THE LARGE MOTION GENERALIZED COORDINATE 68 WPRM ag 
. AND TIME WDOT AND WRESID. * 
x ® 
= ADDITIONAL DEFINITION OF VARIABLES . 
x x 
a HD.......-..-. THE HARTENBERG DENAVIT TRANSFORMATION MATRIX * 
= HDPRM.........-DERIVATIVE OF HD(1) WRT 86 " 
: We wee eeeee eee e+ TRANSFORMATION MATRIX - 
sd Wis ce 0 5s ses 6 »~ TRANSPOSE OF THE TRANSFORMATION MATRIX * 
. WERM.. 2s. 2ae DERIVATIVE OF W WRT 6 i 
= WPRMT......... TRANSPOSE OF THE DERIVATIVE OF W WRT 6 x 
* WOOT... cc ene DERIVATIVE OF W WRT TIME ai 
‘a WOOTT..... .... TRANSPOSE OF THE DERIVATIVE OF W WRT TIME e 
WRESID........ RESIDUAL ACCELERATION TRANSFORMATION > 
* WROUDT = 9 ss ce TRANSPOSE OF THE RESIDUAL ACCELERATION = 
: TRANSFORMATION . 
* | Sy ie Va aoacmra ser VARIOUS KERNELS FOR THE BOOM VOLUME INTEGRAL * 
x x 
® ® 


RRAAAKRKRAARRARARARRARRARARARKRAARARAARARAAARARRAKRAARAARRAARAKRARAKRKKKARKKKHR 


INTEGER TERM 

REAL*8 @(4,4),HD(2,4,4),8PRM(4,4),WDOT(4,4),WRESID(4,4), 
Q,QDOT,Q2D0T, 
WT(4,4),WPRMT(4,4),WDOTT(4,4),WRSIDT(4,4), 
HDPRM(4,4),HD2(4,4),PRD(4,4), 
WD1(4,4),48D2(4,4),4¥D3(4,4),WD4(4,4),8D6(4,4),WD9(4,4), 
WD10(4,4),WD139(4,4),8014(4,4) 


+e ett 
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10 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


DO 


DO 


/BLK15/ Q,QDOT,Q2D0T 

/BLOCK2/ HD 

/BLOCK3/ W 

/BLOCK4/ #D1,WD2,0WD3,WD4,WD6,WD9,WD10,WD13,WD14 
/BLOCK6/ WPRM,WDOT,WRESID 


5 K=1,4 

DO 5 L=1,4 
HDPRM(K,L)=0.0 
HD2(K,L) =HD(2,K,L) 
WDOT(K,L)=0.0 
WRESID(K,L)=0.0 


10 K = 1,4 
HDPRM(2,K)=-HD(1,3,K) 
HDPRM(3,K)= HD(1,2,K) 


CALL MATMUL(PRD,HDPRM,4,4,HD2,4,4,WPRM) 


* THE FOLLOWING MATRIX HAS A QDOT COEFFICIENT 
DO 15 [=1,4 


15 


20 


DO 15 J=1,4 
IF(ABS(WPRM(I,J)).LE.1.E-30) THEN 
wOOT(I,J)=0.D0 
ELSE 
WOOT(I,J)=WPRM(I,J)*QDOT 
ENDIF 


CONTINUE 
* CREATING THE RESIDUAL ACCELERATION TRANSFORMATION MATRIX 
* THIS MATRIX HAS A QDOT**2 COEFFICIENT 


DO 


CONTINUE 


CALL 
CALL 
CALL 
CALL 


CALL 
CALL 
CALL 
CALL 
CALL 


CALL 
CALL 


CALL 
CALL 


END 


20 K = 1,4 
IF(ABS(WDOT(2,K)).LE.1.E-30) THEN 
WRESID(2,K)=0.D0 
ELSE 
WRESID(2,K) =-WDOT(3,K) *QDOT 
ENDIF 
IF(ABS(WDOT(3,K)).LE.1.E-30) THEN 
WRESID(3,K)=0.D0 
ELSE 
WRESID(3,K)=WDOT(2,K) *QDOT 
ENDIF 


TRANSP( TERM, W,4,4,WT) 

TRANSP( TERM, WPRM, 4,4, WPRMT) 
TRANSP( TERM, WDOT,4,4,WDOTT) 
TRANSP( TERM, WRESID, 4,4, WRSIDT) 


MATMUL(PRD,WT,4,4,W,4,4,8D1) 
MATMUL( PRD, WPRMT,4,4,W,4,4,WD2) 
MATMUL(PRD,WDOTT,4,4,W8,4,4,WD3) 
MATMUL(PRD,WRSIDT,4,4,0,4,4,WD4) 
MATMUL( PRD, WPRMT,4,4,WPRM,4,4,W8D6) 


MATMUL(PRD,WT,4,4,WDOT,4,4,WD9) 
MATMUL( PRD, WPRMT,4,4,WDOT,4,4,WD10) 


MATMUL(PRD,WT,4,4,WRESID,4,4,WD13) 
MATMUL(PRD,WPRMT,4,4,WRESID,4,4, WD14) 
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SUBROUTINE INTPRD (TERM, INDEX,WD,X,ROW,COL, PROD) 


RARER RARKRARARAAARRAKKAAKRAAAAARRAARKRARAKRARKRKKAKARKRKEAKRAAKRARRAANAKE 


x 
x 
x 
x 


RRAAAKARKAKAAKAKAKAKAAAKRRARRKARAARKARKRRKKKKKKKRARKKKRAKKRKARKAKRRAREAARKERRAAR 


+++ + + 


x 


THIS SUBROUTINE CREATES THE PRODUCT ARRAY OF THE SHAPE FUNCTIONS * 


AND THE COORDINATE TRANSFORMATION MATRICES W(T) 


INTEGER INDEX, ROW,COL, TERM 
REAL*8 PROD(ROW,COL),R(4),P(4,6),¥4D(4,4), 
A(2),D(2), THETA(2),ALPHA(2),AREA,E,G, 
I¥YY,122,MU(2),GRAV, 
P1(4,6),P1T(6,4),GMAT(4,6),G6T(6,4),61(4,6),G1T(6,4), 
WMAT(4,4),C(4,4),RC(4,1), 
PRD(6,6),X 


COMMON /BLOCK1/ A,D, THETA,ALPHA,AREA,E,G,IYY,12Z2,MU,GRAV 


IF (INDEX .EQ. 11) THEN 
CALL VECTR(X,R) 
PROD(1,1)=0D(1,1) + R(2)*WD(2,1) + R(2)*WD(1,2) 
PROD(1,1)=PROD(1,1)+R(2) **2*WO0(2,2) +1IYY*WD(3,3)/(MU(1) *AREA) 
PROD(1,1)=PROD(1,1)+12Z2*WD(4,4)/(MU(1) *AREA) 

END IF 


IF (INDEX .EQ. 12) THEN 
CALL VECTR(X,R) 
CALL COEFP1(X,P1) 
PROD(1,1)=0D(1,2) *P1(2,1)+R(2)*WD(2,2) *P1(2,1) 
PROD(1,2)=8D(1,3)*P1(3,2)+R(2) *WD(2,3)*P1(3,2) 
PROD(1,3)=WD(1,4)*P1(4,3)+R(2) *WD(2,4) *P1(4,3) 
PROD(1,4)=(-122*W0(4,3) + IYY*WD(3,4))/(MU(.1) *AREA) 
PROD (1,5) =WD( 1,4) *P1(4,5)+R(2)*#0D(2,4) *P1(4,5) 
PROD( 1,6) =8D(1,3) *P1(3,6)+R(2) *8D(2,3)*P1(3,6) 

END IF 


IF (INDEX .EQ. 21) THEN 


CALL COEFP1(X,P1) 
CALL VECTR(X,R) 
PROD(1,1)=WD(2,1)*P1(2,1)+R(2) *WD(2,2)*P1(2,1) 
PROD(2,1)=WD(3,1)*P1(3,2)+R(2) *WO( 3,2) *P1(3,2) 
PROD(3,1)=8D(4,1)*P1(4,3)+R(2)*WD(4,2) *P1(4,3) 
PROD( 4,1) =2(-IZZ2*WD(3,4) + LYY*WD(4,3))/(MU(1) *AREA) 
PROD(5,1)=WD(4,1)*P1(4,5)+R(2) *WD( 4,2) *P1(4,5) 

: PROD(6,1)=24¥D(3,1)*P1(3,6)+R(2)*WD( 3,2) *P1(3,6) 

END IF 


IF (INDEX .EQ. 22) THEN 
CALL COEFP1(X,P1) 
PROD(1,1)=P1(2,1) * WD(2,2) * P1(2,1) 
PROD(1,2)=P1(2,1) * WD(2,3) * P1(3,2) 
PROD(1,3)=P1(2,1) * WD(2,4) * P1(4,3) - 
PROD(1,4) =0.D0 


PROD(1,5)=P1(2,1) * WD(2,4) P1(4,5) 
PROD(1,6)=P1(2,1) * WD(2,3) P1(3,6) 
PROD(2,1)=PROD(1, 2) 

PROD(2,2)=P1(3,2) * WD(3,3) P1(3,2) 
PROD(2,3)=P1(3,2) * WD(3,4) P1(4,3) 
PROD(2,4)=0.D0 

PROD(2,5)=P1(3,2) * WD(3,4) P1(4,5) 
PROD(2,6)=P1(3,2) * WO0(3,3) P1(3,6) 
PROD(3,1)=PROD(1, 3) 

PROD( 3,2) =PROD( 2, 3) 

PROD(3,3)=P1(4,3) * WD(4,4) P1(4,3) 
PROD( 3,4) =0.D0 

PROD(3,5)=P1(4,3) * WD(4,4) P1(4,5) 
PROD(3,6)=P1(4,3) * WD(4,3) P1(3,6) 


PROD(4,1)=0.D0 
PROD(4,2)=0.D0 
PROD (4,3)=0.D0 
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x 


PROD(4,4)=(I2Z * WD(3,3) + IYY * WD(4,3))/(MU(1)*AREA) 


PROD(4,5)=20.D0 
PROD(4,6)=0.D0 
PROD(5,1)=PROD(1,5) 
PROD(5,2)=PROD( 2,5) 
PROD(5, 3) =PROD(3,5) 
PROD(5,4)=0.D0 


PROD(5,5)=P1(4,5) * WD(4,4) * P1(4,5) 
PROD(5,6)=P1(4,5) * WD(4,3) * P1(3,6) 


PROD(6,1)=PROD(1,6) 
PROD(6,2)=PROD( 2.6) 
PROD(6,3)=PROD(3,6) 
PROD(6,4)=0.D0 

PROD(6,5)=PROD(5,6) 


PROD(6,6)=P1(3,6) * WD(3,3) * P1(3,6) 


END IF 


IF (INDEX .EQ. 


CALL MATMUL(PRD,C,4,4,G61,4,6,G1) 
CALL MATMUL(PRD,G1T,6,4,G1, 4,6, PROD) 


END IF 


IF (INDEX .EQ. 


IF (INDEX .EQ. 


1) THEN 
CALL COEFG1(X,G1) 

CALL TRANSP(TERM,G1,4,6,G1T) 
CALL STIFNS(C) 


10) THEN 
CALL VECTR(X,R) 

CALL TRANSP(TERM,R,4,1,PROD) 
END IF 


20) THEN 


CALL COEFP1(X,P1) 


CALL TRANSP (TERM, P1,4,6.PROD) 


END IF 
END 


SUBROUTINE COEFG1 (X,G1) 


e 
RREARRAARRAAKAKAARAARRRRAAARRARRARRARRKAARARARARRARRARARARARRARRAARAARARARAARKRARARRAARARR 


x 
x 
x 


THIS SUBROUTINE CALCULATES THE VALUES FOR THE GAMA1 COEFFICIENT 


x 
% 
2 


RAAARARARAARAARAARARAARARARRARRARARARARRARARARARRARARARARARRARARARARARARARARRAARARAAEKRE 


10 


+ 
+ 


REAL*8 A(2),0(2),THETA(2),ALPHA(2),AREA,E,G, 
IYY,I122,MU(2),GRAV, 


G1(4,6),L,X,K 


COMMON /BLOCK1/ A,D, THETA, ALPHA, AREA,E,G,IYY,122,4U,GRAV 


L=A(2) 
K = X/L 


DO 10 J=1.4 


DO 10 M=1, 


6 


G1(J,M)=0.D0 


G1(1,1) 

G1(2, 3) 

G1(2,5) 

G1(3,2) 

G1(3,6) 

G1(4,4) 
END 


ae 


=(12.*K + °6.)/0222. 


(6.*K + 4.)/L 


G1(2,3) 
G1(2,5) 
G1(1,1) 
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SUBROUTINE STIFNS (5S) 


RRMA RARARARRARAARAAARRARARARARARAARARARAERARARAAARARAARAERAEEK 


x x 


if THIS SUBROUTINE CREATES THE STIFFNESS MATRICES FOR EACH LINK . 
. GIVEN THE MATERIAL AND AREA DATA FROM COMMON BLOCK7 : 
x x 
PrrerrrrrTerrrrrrrrerrrrrrrrrrccrrrrrcrceerCCeOCeecCecrcceroe cere re eee rere errr se. 
REAL*8 A(2),D(2),THETA(2),ALPHA(2),AREA,E.G, 
+ ITYY,I12Z2,MU(2),GRAV,S(4,4),P 


COMMON /BLOCK1/ A,0O, THETA, ALPHA, AREA,E,G,IYY,122,MU,GRAV 


DO 10 I=1,4 

DO 10 J=1,4 
10 S(I,J)=0.00 

P = ITYY+I22Z 
S(1,1) = E*AREA 
S(2,2) = E*IyYy 
S$(3,3) = E*I2Z2 
S(4,4) = G*P 

END 


SUBROUTINE VECTR (X,R) 


RARAEARARARAARAKRARRAAARARARARAARAEARARARAAARAARARARARARARRARARARAARARARARARRARRRARRARRAR 


2 ® 
: THIS SUBROUTINE CALCULATES THE POSITION VECTOR R : 
® s 


RRRAAAAAARARAARAAAAAAAARAARAAARARAARRARRAEAARAARARAARARARARRAARARAARAAARAARARRRAARAARARRAR 


REAL*8 A(2),0(2),THETA(2),ALPHA(2),AREA,E,G, 
+ TYY,12Z2,MU(2),GRAV, 
+ R(4,1),L,X 


COMMON /BLOCK1/ A,O, THETA, ALPHA, AREA,E,G,1IYY,122,4U,GRAV 
L=A(2) 
* CREATING THE LOCAL R POSITION VECTOR 


DO 10 J=1,4 
10 R(J,1)=0.00 : 


R(1,1) = 1.00 
R(2,1) = X 


END 


SUBROUTINE COEFP1 (X,P1). 


RAAARAARARRAAARARRARRARRARARAARARARAARARARAARRAARAARARAARAARARARAAARARAARRAARRARARRRAARNR 
R 2 
* THIS SUBROUTINE WILL CALCULATE THE VALUES FOR THE PHI1 COEFFICIENT * 
2 2 
RARRRAAARRAARAAARARARRARARARRAAARARARARARAARARAARRARAAARRAARRAAARARARRAAKRKARKAKRRAR 


REAL*8 A(2),0(2),THETA(2),ALPHA(2),AREA,E,G, 

+ TYY,I122,MU(2),GRAV, 

+ P1(4,6),L,X,K 

COMMON /BLOCK1/ A,0,THETA,ALPHA, AREA,E,G,IYY,IZZ,MU,GRAV 


L=A(2) 
K = X/L 
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DO 10 J=1,4 
DO 10 M=1,6 


10 P1(J,M)=0.D0 
Pl(2,1) = 1. +K 
P1l(3,2) 's= 2. "°K**O = 9.7 N oe cee 
P1(4,3) = P1(3,2) 
P1(4,5) ‘= L*K*(K**2 + 2.*K + 71.) 
P1(3,6) = P1(4,5) 


END 


SUBROUTINE BINT (TERM,XL,XU,Z,Y,NR,NC, INDEX,DD, MAT) 
RRAARARAAAAARARARARAAARRARAAAARAARAAREAARAKRAARARRARARKRRAEAKRAKRARARAARAKRAHAKAAKRARKRKRARKR 
. GIVEN THE TEMPORARY MATRICES Z, AND Y THE INDEX # 

: CORRESPONDING TO INTEGRAL DEFINITIONS, THE DIMENSIONS OF THE 

. RESULTANT MATRIX PRODUCT (ROW,COL), AND THE W FUNCTION INDEX (D) 
* THIS SUBROUTINE WILL INTEGRATE THE RESULTANT PRODUCT MATRIX 

: OVER THE LENGTH OF THE LINK USING GAUSSIAN QUADRATURE METHODS 

x AN ASSUMPTION OF CONSTANT WEIGHT DENSITY AND CROSS SECTIONAL 

. AREA HAS BEEN MADE WRT THE LINK/ACTUATOR VOLUME INTEGRAL. 

2 


WRRRARAAKAAAAAAARAEAAAAAARARRAARAARAAARARAAAAAAARAARARKAAARARKRAARAAARARAAKAAK 


» » © &® © © DH 


INTEGER NR,NC, INDEX, TERM 
REAL*8 A(2),D(2),THETA(2),ALPHA(2),AREA,E,G, 

+ IYY,IZ2,MU(2),GRAV, 

+ AA,B,C,XL, XU,2(NR,NC), ¥(NR,NC),MAT(NR,NC),DD(4,4) 
COMMON /BLOCK1/ A.D, THETA, ALPHA, AREA,E,G,1IYY,122,MU.GRAV 
AA=.5D0*( XU+XL) 

B=XU~XL 
C=. 48014492824876812D0*B 
CALL INTPRD (TERM, INDEX,D0D, AA+C.NR,NC, Y) 
CALL INTPRD (TERM, INDEX,DD,AA-C,NR,NC,2) 
DO 10 I=1,NR 
DO 10 J=1,NC 
10 MAT(I,J)=.50614268145188130D-1*(Y¥(1I,J)+Z(1,J)) 

C=. 39833323870681337D0*B 
CALL INTPRD (TERM, INDEX,DD,AA+C,NR,NC,Y) 

CALL INTPRD (TERM, INDEX,DD,AA-C,NR,NC,2Z) 

DO 20 I=1,NR 
DO 20 J=1,NC 
20 MAT(I,J)=MAT(I,J)+.11119051722668724D0*(Y(I,J)+2(1,J)) 
C=. 262766 20495816449D0*B 
CALL INTPRD (TERM, INDEX,DD,AA+C,NR,NC, Y) 
CALL INTPRD (TERM, INDEX,DD,AA-C,NR,NC,2) 
DO 30 I=1,NR 
DO 30 J=1,NC 
30 MAT(I,J)=MAT(I,J)+. 15685332293894364D0*(Y¥(I,J)+Z2(1I,J)) 
C=.9171732124782490D~-1"*B 
CALL INTPRD (TERM, INDEX,DD,AA+C,NR.NC, Y) 
CALL INTPRD (TERM, INDEX,DD, AA-C,NR,NC,2Z) 
DO 40 I=1,NR 
DO 40 J=1,NC 
40 MAT(I,J)=B*(MAT(I,J)+.18134189168918099D0*(¥(I,J)+Z(1I,J))) 
DO 50 I=1,NR 
DO 50 J=1,NC 
IF ( INDEX. EQ. 1) THEN 
MAT(I,J)= MAT(I,J) 


ELSE 
MAT(I,J)=(MU(1)*AREA) * MAT(I,J) 
ENDIF 
50 CONTINUE 
RETURN 
END 
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x 
x 
x 
x 


SUBROUTINE DFORM(M,K,UTERM, DC) 


RRAAKAAARAARARARARRARRAARRARAARKARARARARARKRRARRARRRARAAKRAAAKRARAARARARARKARRARKRAARREER 


THIS PROGRAM CALCULATES THE DEFORMATION TRANSFORMATION MATRIX 
FOR THE KINETIC AND POTENTIAL ENERGY CONTRIBUTION OF POINT LOAD 


x 
x 
x 
® 


RARAARARAAARARKRARARKRAAKKRARAKRAAAARRKRARAAKARAKAAKAAKAARAKAARRARARARARRARAARARRRARARK 


10 


REAL*8 UTERM(6),DC(4,4 


DO 10 I=1,4 
DO 10 J=1,4 
DC(I,J)=0.00 


IF(M.EQ.1) THEN 
DC(1,1)=1.00 
DC(2,2)=1.D0 
0C(3,3)=1.00 
DC(4,4)=1.D0 

ENDIF 


TFUCK- EQ. 1). .OR. (K.EQ 
O0C(2,1)=UTERM(1) 
ENDIF 


TFRCCK7EQ-2)' -OR. (K.EQ 
DC(3,1)=UTERM( 2) 
ENDIF 


IF((K.EQ.3) .OR. (K.EQ 
DC( 4,1) =UTERM(3) 
ENDIF 


~IF((K.EQ.4) .OR. (K.EQ 


DC(3,4)=-UTERM(4) 
DC(4,3)=UTERM( 4) 
ENDIF 
IF((K.EQ.5) .OR. (K.EQ 
O0C(2,4) =-UTERM(5) 
0C(4,2)= UTERM(5) 
ENDIF 


IF((K.EQ.6) .OR. (K.EQ 
DC(2,3)=-UTERM(6) 
DC(3,2)=UTERM(6) 

ENDIF 

RETURN 

END 


SUBROUTINE CCOEF(EPSLON,IOTA,H) 


) 


-7)) THEN 


-7)) THEN 


-7)) THEN 


-7)) THEN 


-7)) THEN 


-7)) THEN 


RRARARARAARAARARAARARRARARARRARARARARARRARAAARARARAARARAAKRARARKRAKRAKAAAAKKRAAKAKAAKRAETN 


x 
x 
® 


FORMULATING THE COEFFICIENTS USED IN THE SEQUENTIAL INTEGRATION 
ADDITIONAL DEFINITION OF VARIABLES 
Po) casei sree THE COEFFICIENTS FOR NUMERICAL INTEGRATION 


R 
R 


x 


RRARAKRARAAARARAARARARARRARRARAKRARARARARARAAKRAARARRARARARAAARKARRAKRARARAARRARKARRER 


+ 


REAL*8 EPSLON,IOTA,H, 
A0,A1,A2,A3,A4, 


COMMON /BLOCK8/ AO,Al1,A2,A3,A4,A5,A6,A7 


AO = 1.D0/EPSLON/(H*H) 
Al =IOTA/EPSLON/H 


A5,A6,A7 


DO) 


A2 = 1.D0/EPSLON/H 

A3 = 0.5/EPSLON - 1.D0 

A4 = IOTA/EPSLON - 1.D0 
AS = H/2.D0 * (A4 - 1. 

A6 = H*(1.D0 - IOTA) 

A7 = IOTA * H 

END 
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SUBROUTINE TRANSP (TERM,M1,R,C,M1T) 


RARARARAAARARAAARAARARAARARAARRARRARRRARRAAERARAAARARRRRRRARARRRERERRRERA 
RX 


2 


2 


GIVEN A MATRIX M THIS SUBROUTINE WILL RETURN THE TRANSPOS * 
wl OF THE MATRIX, MT aes i 
® 


x 
MARAARARAAAARARAAKKAARRARAKARAREAA AAA RRR RARER ERERE 


INTEGER R,C, TERM 
REAL*8 M1(R,C), M1T(C,R) 


DO 10 J=1,C 
DO 10 K=1,R 
M1T(J,K)=M1(K,J) 
10 CONTINUE 
END 


SUBROUTINE RINT (1I,4,NOME,J,K) 


RARARARARRARRAARAARARKRRARAAARKRARAARARAAARAAARAARRAARAARAARARAARARARAAAARARARAARKRARRRRAE 


a : a 
* THIS SUBROUTINE WILL PRINT ANY OF THE ARRAYS USED, IT IS * 
” PRIMARILY USED FOR CHECKING THE RETURNED VALUES. i! 
* ” 


RARRARAAAAAAEARAARARREAAARAARAAAARARRAARARARAKRRARRARARARARAARARRAAAKRAARAARRER 


INTEGER I,J,K,L,N,0 
REAL*8 M(I,J,K) 
CHARACTER*6 NAME, NOME 


NAME=! : 
NAME=NOME 
DO 10 L #2 1,I 
WRITE(6,1000) ‘ARRAY ‘,NAME,' FOR I = ',L 
DO 10 N = l,J 
WRITE(6,1100) (M(L,N,0O),0#1,X) 
10 CONTINUE 
1000 FORMAT (1X,A8,:A6,A9,1I1) 
1100 FORMAT (1X,7(: 2X,E9.3)) 
* END 


SUBROUTINE REFORM 


RRARRARRARAAARAKRAAARARARAARAARAARRARRRAAEARARRARAAAARRARAARRAAARAAARRAEAKRAERARAREEAR 


x x 
i REFORMULATING THE COEFFICIENT MATRICES OF THE SMALL MOTION = 

m EQUATION OF MOTION FOR SOLUTION BY SEQUENTIAL INTEGRATION . 
x R 


RRAARARARAARRRARAAERARERARARAERAEAARARRAEAAAAAARARARRAKRAARAAEAARAARARAARARAARKRAARARRARARRARRAR 


REAL*8 MQQ,.MQN(6),MNN(6,6),GN(6,6),KN(6,6),FN(6),FQ, 
+ MN(6,6),FNN(6),FQTEMP,MQNTMP(6),TEMP1(6,6),PRD(6,6) 


COMMON /BLOCK5/ MQQ,MQN.MNN,GN,KN,FN, FQ 
COMMON /BLOCK7/ MN,FNN 


FQTEMP=(1.D0/MQQ) *FQ 
DO 10 [=1,6 
10 FNN(I)=FN(I)-MQN(I) *FQTEMP 
DO 20 I=1,6 
MONTMP(I)=(1.D0/MQQ) *MQN(TI) 
20 CONTINUE 
CALL MATMUL(PRD,MQN,6,1,MQNTMP,1,6, TEMP1) 
DO 30 I[=1,6 
DO 30 J=1,6 
30 MN(I,J)=MNN(I,J)-TEMP1(I,J) 


END 
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SUBROUTINE TRANSF( TIME) 


RARRAKAHHAAARRRAEERRRAAARRAHRARHRAARTRARARARARERRARARRRRERRRERAERRRRRRAEE 


BODY FIXED COORDINATE SYSTEMS 


#4 e he 


COMMON /BLK15/ Q,QDOT,Q2D0T 
COMMON /BLOCK2/ HD 


* CREATING THE HARTENBERG-DENAVIT MATRIX FOR EACH I 


DO 10 I[=1,4 
DO 10 J=1,4 
HD(1,I,J)=0.D0 
10 HD(2,1I,J)=0.D0 
HD(1,1,1) = 1.00 
HD(1,2,2) = CoS(Q) 
HD(1,2,4) = SIN(Q) 
HD(1,3,2) = SIN(Q) 
HD(1,3,4) = -COS(Q) 
HD(1,4,3) = 1.D0 
HD(2,1,1) = 1. 
HD(2,2,1) = A(2)*COS(THETA(2) ) 
HD(2,2,2) = COS(THETA(2)) 
HD(2,2,3) = -SIN(THETA(2) ) 
HD(2.3,1) = A(2)*SIN(THETA(2) ) 
HD(2,3,2) = SIN(THETA(2)) 
HD(2,3,3) = COS(THETA(2)) 
HD(2,4,4) = 1.D0 
RETURN 
END 


SUBROUTINE MATMUL (PRD,MAT1,R1,C1,MAT2,R2,C2,PROD) 


THIS SUBROUTINE CREATES THE HARTENBERG-DENAVIT MATRIX FOR EACH 
JOINT AND THE TRANSFORMATION MATRICES RELATING EACH OF THE JOINT 


ADDITIONAL DEFINITION OF VARIABLES 
HD recess ».++.THE HARTENBERG DENAVIT TRANSFORMATION MATRIX 
RRARAKRAAAAARARRARRARAREEEEAR ARERR AREER EREREAR ARERR RRR ED 
REAL*8 A(2),0(2),THETA(2),ALPHA(2),AREA,.E,G, 
a IYY,122,MU(2),GRAV,HD(2,4,4),Q,QDO0T,Q2D0T, TIME 
COMMON /BLOCK1/ A,D,THETA,ALPHA, AREA,E,.G,IYY,1IZ22,MU,GRAV 


x 
x 
® 
x 
® 
x 


RRRAARARAARARARARARAARAAARRARARARARAARAAARARARARARARRARARARARARAAARKRAARRKARARARRARARKRKRRKRARARAEKRRHRRR 


® 


ad THIS PROGRAM RECEIVES TWO MATRICES TO BE MULTIPLIED TOGETHER 
x ENSURES PROPER DIMENSIONALITY AND RETURNS THE PRODUCT MATRIX 


® 


r 
x 
x 
2 


RRAAAAAHKRKAARRAKRAKRRARRARKRAEARAARAAAARRARAAARAAARARAKARAAAAAAAKARARRARKKRKRAKRKRAR 


INTEGER R1,C1,R2,C2 


REAL*8 MAT1(R1,C1),MAT2(R2,C2),PROD(R1,C2),PRD(R1,C2),SUM 
* CHECKING THE DIMENSIONALITY OF THE MATRICES TO BE MULTIPLIED 
IF (C1.NE.R2) PRINT*, "ERROR IN MULTIPLICTION DIMENSIONS’ 


* MULTIPLYING THE MATRICES TOGETHER 


DO 30 J=1,R1 
DO 20 K=1,C2 
SUM=0.0 
DO 10 M=1,Cl 


IF(ABS(MAT1(J,M)).LE.1.E-30) MAT1(J,M)=0.D0 


IF(ABS(MAT2(M,K)).LE.1.E-30) MAT2(M,K)=0.D0 


CONTINUE 
CONTINUE 
SUM = SUM + MAT1(J,4M) 
10 CONTINUE 
PRD(J,K)=SUM 
20 CONTINUE 


30 CONTINUE 
DO 50 I=1,R1 
DO 50 J=1,C2 
50 PROD(I,J)=PRDO(I,J) 


END 
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* MAT2(M,X) 


SUBROUTINE WHATRX 


RRARRRARARARARAARARARAARARARARAARARRARARARAARRARARARARARARARARARARRAARARARARARAARARARARARKRAAARAR 


2 ® 
* THIS SUBROUTINE MULTIPLIES THE HARTENBERG-DENAVIT MATRICES ‘ 
* TOGETHER TO OBTAIN THE COORDINATE TRANSFORMATION MATRIX * 
2 ® 
* ADDITIONAL DEFINITION OF VARIABLES * 
Rx n 
* (ar re THE HARTENBERG DENAVIT TRANSFORMATION MATRIX * 
* vee weeeeeess TRANSFORMATION MATRIX * 
R 2 
2 


RRARAARRARARRARARARARAARARAARAARARAARAAARAAAARARARARAAAAARATAAARAARARRARARARRAARAARARARARAARRAR 


REAL*S W(4,4),HD(2,4,4),WSUM(4,4),5UM 


COMMON /BLOCK2/ HD 
COMMON /BLOCK3/ WwW 


DO 10 L=1,4 
DO 10 M=1,4 
WSUM(L,M)=0.D0 
10 W(L,M)=0.D0 
DO 40 J=1,4 
DO 30 K=1,4 
SUM=0.0 


DO 20 M=1,4 
SUM = SUM + HD(1,J,4) * HD(2,M,X) 


20 CONTINUE 
WSUM(J,K)=SUM 
30 CONTINUE 
40 CONTINUE 
DO 50 [=1,4 
DO 50 J=1,4 
50 W(I,J) =WSUM(I,J) 
END : 
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Initialize 
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time step, and 
finish time 
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call CCOEF 
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time depend 
integration 
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call INTSML 
solves for UV 
call INTLRG 
solves for 8 


print output 




















call TSML 
implicit 
call SMLCOF 
Mnn,Gn,Kn,Fn 
call JINTLRG 
explicit. 
call LRGCOF 
Mqq,Mqn,Fq 
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Figure 15 Simulation Code Flow Chart 
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